Model selection for Ch 1 - marsh data

*NOTES:*

  1. Decided we should NOT include site as a random effect for all beach seine models because the low # of sites and it performed worse than regular GLMs for all.
  1. Early on in model selection we determined that our data has a negative binomial distribution. All models started with Poisson distributions and all were overdispersed, which was corrected by the neg binomial distribution. All models here use negative binomial distribution with a log link.
  1. We used AICc to rank all models in dredge, as it penalizes the strength of likelihood more for very small sample sizes, and is more adaptive to moderate sample sizes, so that we could apply the same information criterion across all models (Harrison et al 2018; Brewer et al 2016) – could explain better…
  1. We used model averaging following a delta AICc of less than 4

Chinook beach seine model

response = average abundance of Chinook salmon [per beach seine site]
This is 1 of 4 separate marsh models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from beach seine data

# load Beach seine data aggregated by species/set; note this contains
# unidentified larval fish
b <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/beach.catch.csv")

summary(b)
##       Year               Date       Temp.surf        DO.surf      
##  Min.   :0.0000   5/17/2017: 42   Min.   : 3.43   Min.   : 48.80  
##  1st Qu.:0.0000   30-May-16: 35   1st Qu.:11.68   1st Qu.: 90.90  
##  Median :0.0000   5/31/2017: 32   Median :13.30   Median :103.20  
##  Mean   :0.4817   6/16/2017: 31   Mean   :13.36   Mean   : 99.55  
##  3rd Qu.:1.0000   23-Sep-16: 24   3rd Qu.:16.30   3rd Qu.:112.80  
##  Max.   :1.0000   7/12/2017: 24   Max.   :18.32   Max.   :136.60  
##                   (Other)  :632                                   
##    DOmg.surf        pH.surf          Sal.surf         Temp.bot    
##  Min.   : 4.68   Min.   : 6.440   Min.   : 0.000   Min.   : 4.15  
##  1st Qu.: 9.07   1st Qu.: 7.710   1st Qu.: 0.110   1st Qu.:11.58  
##  Median :10.74   Median : 7.905   Median : 0.890   Median :13.20  
##  Mean   :10.46   Mean   : 8.059   Mean   : 2.326   Mean   :13.36  
##  3rd Qu.:11.72   3rd Qu.: 8.150   3rd Qu.: 3.270   3rd Qu.:16.30  
##  Max.   :15.50   Max.   :10.460   Max.   :19.400   Max.   :18.32  
##                                                                   
##      DO.bot          DOmg.bot         pH.bot         Sal.bot      
##  Min.   : 48.50   Min.   : 4.68   Min.   :6.440   Min.   : 0.050  
##  1st Qu.: 93.20   1st Qu.: 9.21   1st Qu.:7.750   1st Qu.: 0.120  
##  Median :101.40   Median :10.49   Median :7.895   Median : 0.850  
##  Mean   : 99.13   Mean   :10.40   Mean   :8.021   Mean   : 2.416  
##  3rd Qu.:112.50   3rd Qu.:11.54   3rd Qu.:8.150   3rd Qu.: 3.340  
##  Max.   :140.30   Max.   :15.50   Max.   :9.890   Max.   :20.750  
##                                                                   
##      J.date         SIN_TIME          COS_TIME           Round       
##  Min.   : 64.0   Min.   :-0.9929   Min.   :-0.9990   Min.   : 1.000  
##  1st Qu.:124.0   1st Qu.:-0.1436   1st Qu.:-0.9641   1st Qu.: 4.000  
##  Median :151.0   Median : 0.5176   Median :-0.7075   Median : 6.000  
##  Mean   :157.2   Mean   : 0.3478   Mean   :-0.6395   Mean   : 5.901  
##  3rd Qu.:191.0   3rd Qu.: 0.8460   3rd Qu.:-0.4431   3rd Qu.: 8.000  
##  Max.   :282.0   Max.   : 1.0000   Max.   : 0.4527   Max.   :10.000  
##                                                                      
##   Habitat    Site          Set                            Species   
##  Marsh:820   M1:155   Min.   :1.000   Three-spined stickleback:161  
##              M2:179   1st Qu.:1.000   Chinook                 :124  
##              M3:171   Median :2.000   Peamouth chub           :119  
##              M4:172   Mean   :2.059   Starry flounder         : 98  
##              M5:143   3rd Qu.:3.000   Staghorn sculpin        : 74  
##                       Max.   :7.000   Prickly sculpin         : 50  
##                                       (Other)                 :194  
##     Life_stage    abundance            class         round       
##  0       : 28   Min.   :  0.00   0        : 28   Min.   : 0.500  
##  adult   :176   1st Qu.:  1.00   migratory:201   1st Qu.: 4.000  
##  juvenile:474   Median :  2.00   resident :588   Median : 6.000  
##  NA's    :142   Mean   :  8.69   NA's     :  3   Mean   : 5.757  
##                 3rd Qu.:  6.00                   3rd Qu.: 8.000  
##                 Max.   :509.00                   Max.   :10.000  
##                                                                  
##        month    
##  May      :210  
##  June     :163  
##  April    :160  
##  July     :119  
##  March    : 81  
##  September: 45  
##  (Other)  : 42
sapply(b, sd)
##       Year       Date  Temp.surf    DO.surf  DOmg.surf    pH.surf 
##  0.4999702 16.0088878  3.2129387 19.5469879  2.0540975  0.6425623 
##   Sal.surf   Temp.bot     DO.bot   DOmg.bot     pH.bot    Sal.bot 
##  3.3653542  3.1011628 19.8481953  2.1137164  0.5439179  3.6501406 
##     J.date   SIN_TIME   COS_TIME      Round    Habitat       Site 
## 47.0456528  0.5773336  0.3706291  2.5508623  0.0000000  1.3720692 
##        Set    Species Life_stage  abundance      class      round 
##  0.9564320  6.3846921         NA 29.7143774         NA  2.2935398 
##      month 
##  2.0433562
### Separate into groups; automatically removes 3 x unidentified larval and 28
### x '0' abundance rows 1
b.1 <- b[which(b$Species == "Chinook"), ]

### standardize continuous variables
library(robustHD)
# standardize variables to be centered on the mean (mean becomes 0) using
# the standardize function from robustHD
b.1$s.temp <- standardize(b.1$Temp.surf, centerFun = mean, scaleFun = sd)
b.1$s.sal <- standardize(b.1$Sal.surf, centerFun = mean, scaleFun = sd)
b.1$s.do <- standardize(b.1$DOmg.surf, centerFun = mean, scaleFun = sd)
b.1$s.pH <- standardize(b.1$pH.surf, centerFun = mean, scaleFun = sd)
b.1$s.J.date <- standardize(b.1$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
b.1$j2 <- b.1$s.J.date^2

# 1 summarize by site-day
library(plyr)
b.chin <- ddply(b.1, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.chin)

Habitat variables

Load habitat variables – these are habitat characteristics measured at each marsh site once. See “habitat characteristics” section of model.b3.html for full descriptions

## Add marsh habitat variables for each site
hab <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/site.char.marsh.csv")
hab2 <- hab[, c(1, 6, 8:10, 13, 19, 20)]

summary(hab2)
##  Site       stnwidth        vegheight        shtdensity       shtdenhigh  
##    :13   Min.   : 18.80   Min.   :0.3235   Min.   : 176.0   Min.   : 752  
##  M1: 4   1st Qu.: 24.30   1st Qu.:0.6747   1st Qu.: 344.0   1st Qu.:1664  
##  M2: 1   Median : 53.90   Median :1.2538   Median : 709.3   Median :1837  
##  M3: 4   Mean   : 51.26   Mean   :1.2451   Mean   : 812.2   Mean   :1843  
##  M4: 4   3rd Qu.: 60.30   3rd Qu.:1.5718   3rd Qu.:1392.0   3rd Qu.:1933  
##  M5: 4   Max.   :170.07   Max.   :2.5837   Max.   :1792.0   Max.   :2768  
##          NA's   :13       NA's   :13       NA's   :13       NA's   :13    
##      tcelev          angbank          meanturb    
##  Min.   :0.3440   Min.   : 4.000   Min.   :29.87  
##  1st Qu.:0.6390   1st Qu.: 8.531   1st Qu.:33.67  
##  Median :0.8470   Median : 9.500   Median :37.47  
##  Mean   :0.8408   Mean   :12.144   Mean   :35.95  
##  3rd Qu.:1.0500   3rd Qu.:12.917   3rd Qu.:41.02  
##  Max.   :1.3590   Max.   :31.000   Max.   :41.52  
##  NA's   :13       NA's   :13       NA's   :13
sapply(hab2, sd, na.rm = TRUE)
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
##   Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
##        Site    stnwidth   vegheight  shtdensity  shtdenhigh      tcelev 
##   1.9546584  37.6011569   0.6223313 552.9042986 475.4831598   0.2944439 
##     angbank    meanturb 
##   6.6741386   4.5283513
b.chin.hab <- b.chin
b.chin.hab$stnwidth <- hab2[match(b.chin.hab$Site, hab$Site), 2]  #add stnwidth to b.chin.hab by matching with Site #. Then repeat with all other habitat variables  
b.chin.hab$vegelev <- hab2[match(b.chin.hab$Site, hab$Site), 3]
b.chin.hab$shtdensity <- hab2[match(b.chin.hab$Site, hab$Site), 4]
b.chin.hab$shtdenhigh <- hab2[match(b.chin.hab$Site, hab$Site), 5]
b.chin.hab$tcelev <- hab2[match(b.chin.hab$Site, hab$Site), 6]
b.chin.hab$angbank <- hab2[match(b.chin.hab$Site, hab$Site), 7]
b.chin.hab$meanturb <- hab2[match(b.chin.hab$Site, hab$Site), 8]

## standardize habitat variables using 'standardize' function from robustHD
## package
b.chin.hab$stnwidth <- standardize(b.chin.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.chin.hab$stnwidth <- as.numeric(b.chin.hab$stnwidth)
b.chin.hab$vegelev <- standardize(b.chin.hab$vegelev, centerFun = mean, scaleFun = sd)
b.chin.hab$vegelev <- as.numeric(b.chin.hab$vegelev)
b.chin.hab$shtdensity <- standardize(b.chin.hab$shtdensity, centerFun = mean, 
    scaleFun = sd)
b.chin.hab$shtdensity <- as.numeric(b.chin.hab$shtdensity)
b.chin.hab$shtdenhigh <- standardize(b.chin.hab$shtdenhigh, centerFun = mean, 
    scaleFun = sd)
b.chin.hab$shtdenhigh <- as.numeric(b.chin.hab$shtdenhigh)
b.chin.hab$tcelev <- standardize(b.chin.hab$tcelev, centerFun = mean, scaleFun = sd)
b.chin.hab$tcelev <- as.numeric(b.chin.hab$tcelev)
b.chin.hab$angbank <- standardize(b.chin.hab$angbank, centerFun = mean, scaleFun = sd)
b.chin.hab$angbank <- as.numeric(b.chin.hab$angbank)
b.chin.hab$meanturb <- standardize(b.chin.hab$meanturb, centerFun = mean, scaleFun = sd)
b.chin.hab$meanturb <- as.numeric(b.chin.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors
Note that despite high collinearity of temp, when we tried to remove it the model would not converge. When we ran just temperature and removed julian day and j2, the AIC and deviance explained (r2) were considerably worse. We determined in this case it would be best to leave it in.

library(car)
library(GGally)

# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.chin.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
##     stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank + 
##     meanturb
## 
## Complete :
##          (Intercept)     s.J.date        j2              Year           
## tcelev                 0               0               0               0
## angbank                0               0               0               0
## meanturb               0               0               0               0
##          s.temp          s.sal           s.do            s.pH           
## tcelev                 0               0               0               0
## angbank                0               0               0               0
## meanturb               0               0               0               0
##          stnwidth        vegelev         shtdenhigh      shtdensity     
## tcelev       17695/31146     16003/35591    -19915/37564 1770469/5798310
## angbank    -58882/110101   258521/352474      -1628/3495   -20031/306857
## meanturb      3877/11296    -48133/68263    -14929/23063      -1997/4121

# We determined hab variables for Chinook should be reduced to stnwidth,
# meanturb, and vegelev based on significance and collinearity when doing
# model exploration (see marsh.chinook.Rmd).
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + meanturb, data = b.chin.hab))
## s.J.date       j2     Year   s.temp    s.sal     s.do     s.pH stnwidth 
## 6.572784 1.992816 1.802451 6.767715 2.158658 2.514453 1.512845 1.920460 
##  vegelev meanturb 
## 1.726731 1.642707


## Pearson Corr with all vars
year <- b.chin.hab$Year
jday <- b.chin.hab$s.J.date
j2 <- b.chin.hab$j2
temp <- b.chin.hab$s.temp
do <- b.chin.hab$s.do
ph <- b.chin.hab$s.pH
sal <- b.chin.hab$s.sal
veg <- b.chin.hab$vegelev
turb <- b.chin.hab$meanturb
width <- b.chin.hab$stnwidth
shootden <- b.chin.hab$shtdensity
shtdenhigh <- b.chin.hab$shtdenhigh
elev <- b.chin.hab$tcelev
angbank <- b.chin.hab$angbank

habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width, 
    shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine Chinook salmon dataset")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chinookmarsh.pdf")

Model selection

Full model

The beach seine models performed better without site as a random effect (e.g. see marsh.chinook.Rmd), and only had a maximum of 5 sites in each data set, so we decided not to include site as a RE in the global models.

library(MASS)

## Chinook: full model with top habitat variables
chin <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
    stnwidth + vegelev + meanturb, data = b.chin.hab, na.action = "na.fail")
summary(chin)  #AIC 441.56, deviance explained = 100*((132.69-64.26)/132.69) = 51.57 % ((note this is same calculation that function below uses for pseudo r-squared))
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp + 
##     s.sal + s.do + s.pH + stnwidth + vegelev + meanturb, data = b.chin.hab, 
##     na.action = "na.fail", init.theta = 1.448876434, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4538  -0.9538  -0.4630   0.4737   1.9198  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.857958   0.231350   8.031 9.67e-16 ***
## s.J.date    -1.528136   0.270307  -5.653 1.57e-08 ***
## j2          -0.074091   0.126257  -0.587 0.557320    
## Year         1.043695   0.309056   3.377 0.000733 ***
## s.temp       0.927479   0.274609   3.377 0.000732 ***
## s.sal       -0.097340   0.161957  -0.601 0.547825    
## s.do        -0.115719   0.174303  -0.664 0.506757    
## s.pH         0.210129   0.146060   1.439 0.150249    
## stnwidth     0.372425   0.157667   2.362 0.018172 *  
## vegelev     -0.044583   0.150723  -0.296 0.767389    
## meanturb     0.002958   0.148951   0.020 0.984157    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.4489) family taken to be 1)
## 
##     Null deviance: 132.69  on 62  degrees of freedom
## Residual deviance:  64.26  on 52  degrees of freedom
## AIC: 441.56
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.449 
##           Std. Err.:  0.279 
## 
##  2 x log-likelihood:  -417.559
chin2 <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
    shtdensity + angbank + tcelev, data = b.chin.hab, na.action = "na.fail")
summary(chin2)  #AIC 441.55
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp + 
##     s.sal + s.do + s.pH + shtdensity + angbank + tcelev, data = b.chin.hab, 
##     na.action = "na.fail", init.theta = 1.452371266, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4044  -0.9797  -0.4175   0.5125   2.0049  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.85729    0.23027   8.066 7.28e-16 ***
## s.J.date    -1.51924    0.26942  -5.639 1.71e-08 ***
## j2          -0.07838    0.12493  -0.627 0.530397    
## Year         1.05217    0.31102   3.383 0.000717 ***
## s.temp       0.93394    0.27886   3.349 0.000811 ***
## s.sal       -0.08800    0.16035  -0.549 0.583136    
## s.do        -0.11516    0.17403  -0.662 0.508158    
## s.pH         0.20512    0.14431   1.421 0.155207    
## shtdensity  -0.10413    0.12897  -0.807 0.419426    
## angbank     -0.33802    0.15674  -2.156 0.031046 *  
## tcelev       0.37826    0.16659   2.271 0.023166 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.4524) family taken to be 1)
## 
##     Null deviance: 132.965  on 62  degrees of freedom
## Residual deviance:  64.384  on 52  degrees of freedom
## AIC: 441.55
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.452 
##           Std. Err.:  0.280 
## 
##  2 x log-likelihood:  -417.554
vif(chin2)
##   s.J.date         j2       Year     s.temp      s.sal       s.do 
##   6.155167   1.967313   1.844052   6.557370   2.105972   2.502986 
##       s.pH shtdensity    angbank     tcelev 
##   1.430439   1.251839   1.881877   2.102137
## Now look at the variance inflation factors of this full model
vif(glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
    stnwidth + vegelev + meanturb, data = b.chin.hab))
## s.J.date       j2     Year   s.temp    s.sal     s.do     s.pH stnwidth 
## 6.182463 2.004886 1.816737 6.348249 2.145333 2.501217 1.462386 1.937416 
##  vegelev meanturb 
## 1.706736 1.669399
# Test fit and assumptions of full (global) model - Jday is high (6.18)
library(sjPlot)
## Warning: package 'sjPlot' was built under R version 3.4.3
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.12
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(chin, type = "ma")
## No outliers detected.
## NULL

## --------------------
## Check significance of terms when they entered the model...
## Anova:
## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'

## Analysis of Deviance Table
## 
## Model: Negative Binomial(1.4489), link: log
## 
## Response: abundance
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        62    132.688              
## s.J.date  1   43.330        61     89.359 4.625e-11 ***
## j2        1    1.013        60     88.346 0.3142812    
## Year      1    0.775        59     87.571 0.3786187    
## s.temp    1   10.904        58     76.667 0.0009597 ***
## s.sal     1    0.271        57     76.397 0.6028306    
## s.do      1    1.269        56     75.128 0.2599757    
## s.pH      1    1.343        55     73.784 0.2464322    
## stnwidth  1    9.435        54     64.349 0.0021283 ** 
## vegelev   1    0.089        53     64.260 0.7660493    
## meanturb  1    0.000        52     64.260 0.9849889    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# See pseudo R2 for full model
library(piecewiseSEM)
rsquared(chin, aicc = TRUE)
##    Class            Family Link  n R.squared
## 1 negbin Negative Binomial  log 63  0.515708

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (chin has 63 observations)

library(MuMIn)
library(knitr)
# Generate model set
model.set.chin <- dredge(chin)

# Create model selection table
model_table.chin <- model.sel(model.set.chin)
options(scipen = 7)
names(model_table.chin) <- c("(Int)", "Jday^2", "Mean turb", "DO", "Jday", "pH", 
    "Sal", "Temp", "Width", "Veg. elev", "Year", "df", "LL", "AICc", "delta", 
    "weight")
kable(head(model_table.chin, n = 100), digits = 3)
(Int) Jday^2 Mean turb DO Jday pH Sal Temp Width Veg. elev Year df LL AICc delta weight
729 1.839 NA NA NA -1.449 0.174 NA 0.908 0.323 NA 0.959 7 -209.394 434.824 0.000 0.1139834028
713 1.834 NA NA NA -1.418 NA NA 0.886 0.306 NA 1.003 6 -210.904 435.308 0.484 0.0894893114
969 1.822 NA NA NA -1.441 NA NA 0.905 0.355 -0.130 1.014 7 -210.387 436.810 1.986 0.0422265266
761 1.801 NA NA NA -1.471 0.176 -0.120 0.889 0.335 NA 1.026 8 -209.079 436.825 2.001 0.0419136723
985 1.831 NA NA NA -1.456 0.158 NA 0.913 0.348 -0.069 0.971 8 -209.258 437.183 2.359 0.0350370532
733 1.848 NA NA -0.063 -1.463 0.201 NA 0.895 0.337 NA 0.929 8 -209.283 437.233 2.409 0.0341737131
745 1.800 NA NA NA -1.438 NA -0.111 0.869 0.314 NA 1.064 7 -210.675 437.386 2.562 0.0316545586
715 1.831 NA 0.091 NA -1.396 NA NA 0.866 0.265 NA 1.001 7 -210.683 437.402 2.579 0.0313987993
731 1.838 NA 0.026 NA -1.442 0.169 NA 0.901 0.311 NA 0.961 8 -209.377 437.422 2.598 0.0310983070
730 1.848 -0.016 NA NA -1.456 0.173 NA 0.920 0.323 NA 0.976 8 -209.380 437.427 2.603 0.0310139791
717 1.827 NA NA 0.062 -1.410 NA NA 0.904 0.296 NA 1.026 7 -210.763 437.562 2.738 0.0289933656
714 1.853 -0.034 NA NA -1.432 NA NA 0.912 0.306 NA 1.036 7 -210.845 437.726 2.903 0.0267000550
587 1.884 NA 0.242 NA -1.334 NA NA 0.821 NA NA 0.946 6 -212.644 438.788 3.964 0.0157038322
970 1.842 -0.038 NA NA -1.458 NA NA 0.935 0.356 -0.133 1.053 8 -210.312 439.291 4.467 0.0122144611
973 1.818 NA NA 0.036 -1.435 NA NA 0.914 0.346 -0.121 1.028 8 -210.343 439.352 4.529 0.0118423972
971 1.821 NA 0.041 NA -1.429 NA NA 0.895 0.331 -0.115 1.013 8 -210.350 439.367 4.544 0.0117550191
1001 1.810 NA NA NA -1.446 NA -0.044 0.896 0.353 -0.114 1.037 8 -210.359 439.384 4.560 0.0116560358
765 1.810 NA NA -0.055 -1.482 0.200 -0.116 0.878 0.347 NA 0.996 9 -208.997 439.391 4.567 0.0116169357
762 1.813 -0.026 NA NA -1.484 0.175 -0.125 0.909 0.335 NA 1.057 9 -209.043 439.482 4.658 0.0111019226
763 1.799 NA 0.030 NA -1.463 0.171 -0.121 0.881 0.321 NA 1.028 9 -209.057 439.509 4.686 0.0109485748
747 1.794 NA 0.097 NA -1.417 NA -0.115 0.848 0.272 NA 1.065 8 -210.425 439.516 4.692 0.0109115390
1017 1.803 NA NA NA -1.471 0.171 -0.109 0.892 0.341 -0.020 1.023 9 -209.070 439.537 4.713 0.0108011948
749 1.787 NA NA 0.073 -1.431 NA -0.121 0.889 0.304 NA 1.098 8 -210.481 439.629 4.805 0.0103155662
989 1.840 NA NA -0.069 -1.471 0.187 NA 0.898 0.365 -0.074 0.938 9 -209.126 439.648 4.824 0.0102158399
734 1.886 -0.060 NA -0.109 -1.499 0.217 NA 0.932 0.347 NA 0.971 9 -209.149 439.694 4.871 0.0099816884
746 1.821 -0.043 NA NA -1.459 NA -0.118 0.902 0.316 NA 1.112 8 -210.580 439.827 5.003 0.0093401344
719 1.825 NA 0.081 0.049 -1.392 NA NA 0.882 0.262 NA 1.021 8 -210.597 439.861 5.038 0.0091813130
986 1.843 -0.021 NA NA -1.466 0.157 NA 0.929 0.348 -0.072 0.995 9 -209.234 439.864 5.040 0.0091695253
987 1.831 NA 0.000 NA -1.456 0.159 NA 0.913 0.348 -0.069 0.971 9 -209.258 439.913 5.089 0.0089497627
735 1.846 NA 0.026 -0.063 -1.456 0.196 NA 0.888 0.325 NA 0.930 9 -209.267 439.929 5.106 0.0088745728
716 1.845 -0.025 0.087 NA -1.407 NA NA 0.886 0.267 NA 1.026 8 -210.650 439.968 5.144 0.0087068344
603 1.892 NA 0.214 NA -1.357 0.127 NA 0.836 NA NA 0.919 7 -211.997 440.030 5.206 0.0084387393
585 1.937 NA NA NA -1.373 NA NA 0.855 NA NA 0.904 5 -214.515 440.083 5.259 0.0082185116
732 1.846 -0.015 0.024 NA -1.449 0.169 NA 0.912 0.311 NA 0.976 9 -209.366 440.129 5.305 0.0080338508
718 1.832 -0.008 NA 0.057 -1.414 NA NA 0.908 0.296 NA 1.032 8 -210.761 440.188 5.364 0.0077990801
601 1.941 NA NA NA -1.393 0.163 NA 0.873 NA NA 0.868 6 -213.444 440.388 5.565 0.0070544988
619 1.853 NA 0.247 NA -1.355 NA -0.097 0.807 NA NA 1.001 7 -212.468 440.972 6.148 0.0052690565
591 1.875 NA 0.223 0.070 -1.330 NA NA 0.845 NA NA 0.974 7 -212.483 441.002 6.178 0.0051910483
843 1.883 NA 0.246 NA -1.327 NA NA 0.814 NA 0.039 0.945 7 -212.602 441.241 6.417 0.0046064972
588 1.892 -0.015 0.240 NA -1.341 NA NA 0.834 NA NA 0.960 7 -212.634 441.305 6.481 0.0044615718
589 1.916 NA NA 0.129 -1.359 NA NA 0.893 NA NA 0.955 6 -213.965 441.431 6.607 0.0041892510
217 2.396 NA NA NA -1.017 0.186 NA 0.394 0.289 NA NA 6 -214.040 441.580 6.756 0.0038884300
766 1.850 -0.068 NA -0.106 -1.525 0.218 -0.123 0.919 0.358 NA 1.049 10 -208.824 441.879 7.055 0.0033489426
1002 1.830 -0.042 NA NA -1.467 NA -0.053 0.928 0.353 -0.113 1.085 9 -210.270 441.936 7.112 0.0032541919
1005 1.799 NA NA 0.047 -1.440 NA -0.062 0.904 0.339 -0.095 1.065 9 -210.289 441.974 7.150 0.0031927535
972 1.840 -0.035 0.032 NA -1.447 NA NA 0.924 0.337 -0.121 1.048 9 -210.290 441.976 7.153 0.0031890868
1003 1.803 NA 0.058 NA -1.432 NA -0.066 0.877 0.317 -0.085 1.046 9 -210.292 441.980 7.156 0.0031834838
751 1.784 NA 0.084 0.059 -1.414 NA -0.122 0.868 0.268 NA 1.094 9 -210.300 441.996 7.172 0.0031576250
974 1.837 -0.032 NA 0.014 -1.453 NA NA 0.934 0.352 -0.129 1.052 9 -210.308 442.011 7.188 0.0031338963
975 1.817 NA 0.037 0.033 -1.424 NA NA 0.904 0.325 -0.108 1.026 9 -210.313 442.023 7.199 0.0031158397
990 1.886 -0.075 NA -0.128 -1.518 0.204 NA 0.945 0.383 -0.089 0.992 10 -208.927 442.084 7.260 0.0030217716
748 1.811 -0.034 0.090 NA -1.435 NA -0.120 0.875 0.275 NA 1.102 9 -210.369 442.133 7.310 0.0029482865
767 1.808 NA 0.030 -0.055 -1.474 0.194 -0.117 0.869 0.332 NA 0.999 10 -208.975 442.180 7.357 0.0028797643
1021 1.813 NA NA -0.058 -1.482 0.194 -0.098 0.881 0.356 -0.030 0.990 10 -208.978 442.188 7.364 0.0028693484
859 1.894 NA 0.219 NA -1.348 0.146 NA 0.825 NA 0.087 0.910 8 -211.793 442.254 7.430 0.0027763221
764 1.811 -0.024 0.026 NA -1.476 0.170 -0.125 0.900 0.322 NA 1.057 10 -209.025 442.282 7.458 0.0027378188
635 1.862 NA 0.220 NA -1.377 0.126 -0.096 0.821 NA NA 0.973 8 -211.810 442.286 7.463 0.0027311620
1018 1.815 -0.027 NA NA -1.484 0.169 -0.112 0.913 0.342 -0.022 1.054 10 -209.032 442.294 7.471 0.0027203250
1019 1.800 NA 0.027 NA -1.464 0.169 -0.116 0.883 0.325 -0.008 1.027 10 -209.055 442.341 7.518 0.0026570316
750 1.796 -0.015 NA 0.064 -1.439 NA -0.122 0.898 0.305 NA 1.110 9 -210.473 442.342 7.519 0.0026557726
617 1.914 NA NA NA -1.387 NA -0.076 0.842 NA NA 0.946 6 -214.424 442.347 7.524 0.0026490612
586 1.955 -0.035 NA NA -1.390 NA NA 0.887 NA NA 0.941 6 -214.463 442.426 7.602 0.0025472267
201 2.421 NA NA NA -0.969 NA NA 0.341 0.266 NA NA 5 -215.700 442.453 7.629 0.0025131666
991 1.840 NA -0.002 -0.069 -1.472 0.187 NA 0.899 0.367 -0.075 0.938 10 -209.126 442.482 7.659 0.0024763712
736 1.884 -0.059 0.019 -0.108 -1.493 0.213 NA 0.925 0.338 NA 0.971 10 -209.140 442.511 7.687 0.0024411353
841 1.937 NA NA NA -1.373 NA NA 0.855 NA 0.001 0.904 6 -214.515 442.530 7.707 0.0024174596
720 1.828 -0.005 0.081 0.046 -1.394 NA NA 0.885 0.262 NA 1.024 9 -210.597 442.589 7.766 0.0023472717
604 1.898 -0.011 0.213 NA -1.363 0.127 NA 0.845 NA NA 0.930 8 -211.991 442.650 7.826 0.0022775518
607 1.891 NA 0.212 0.010 -1.356 0.123 NA 0.839 NA NA 0.923 8 -211.994 442.655 7.831 0.0022715015
857 1.943 NA NA NA -1.388 0.179 NA 0.869 NA 0.067 0.862 7 -213.324 442.684 7.860 0.0022391394
988 1.843 -0.022 -0.005 NA -1.467 0.157 NA 0.931 0.351 -0.073 0.995 10 -209.233 442.698 7.874 0.0022234472
633 1.918 NA NA NA -1.407 0.163 -0.077 0.860 NA NA 0.911 7 -213.338 442.712 7.889 0.0022071363
605 1.932 NA NA 0.054 -1.385 0.140 NA 0.885 NA NA 0.894 7 -213.370 442.777 7.953 0.0021375320
602 1.953 -0.024 NA NA -1.406 0.162 NA 0.895 NA NA 0.896 7 -213.419 442.874 8.050 0.0020356294
137 2.438 NA NA NA -0.708 NA NA NA 0.269 NA NA 4 -217.175 443.040 8.216 0.0018740150
221 2.377 NA NA -0.129 -1.076 0.240 NA 0.398 0.326 NA NA 7 -213.558 443.152 8.328 0.0017715904
875 1.837 NA 0.261 NA -1.348 NA -0.148 0.782 NA 0.093 1.027 8 -212.267 443.201 8.377 0.0017288359
623 1.841 NA 0.228 0.078 -1.352 NA -0.104 0.834 NA NA 1.037 8 -212.270 443.207 8.383 0.0017234459
153 2.420 NA NA NA -0.715 0.158 NA NA 0.287 NA NA 5 -216.122 443.298 8.474 0.0016472706
218 2.307 0.077 NA NA -1.021 0.193 NA 0.376 0.292 NA NA 7 -213.726 443.488 8.665 0.0014974741
847 1.874 NA 0.228 0.075 -1.322 NA NA 0.837 NA 0.047 0.974 8 -212.421 443.509 8.685 0.0014822040
620 1.863 -0.019 0.246 NA -1.364 NA -0.098 0.824 NA NA 1.020 8 -212.452 443.570 8.747 0.0014372034
592 1.863 0.019 0.223 0.082 -1.320 NA NA 0.833 NA NA 0.960 8 -212.470 443.606 8.782 0.0014118203
621 1.886 NA NA 0.136 -1.376 NA -0.092 0.881 NA NA 1.010 7 -213.816 443.669 8.845 0.0013682316
891 1.837 NA 0.233 NA -1.378 0.165 -0.187 0.788 NA 0.168 1.007 9 -211.214 443.823 9.000 0.0012664062
844 1.891 -0.012 0.245 NA -1.333 NA NA 0.824 NA 0.038 0.957 8 -212.595 443.857 9.033 0.0012453465
590 1.900 0.025 NA 0.144 -1.344 NA NA 0.874 NA NA 0.934 7 -213.943 443.923 9.099 0.0012050360
845 1.916 NA NA 0.132 -1.356 NA NA 0.891 NA 0.022 0.956 7 -213.952 443.940 9.116 0.0011948853
473 2.395 NA NA NA -1.021 0.176 NA 0.395 0.304 -0.041 NA 7 -213.996 444.028 9.205 0.0011430536
219 2.395 NA 0.014 NA -1.013 0.183 NA 0.389 0.283 NA NA 7 -214.035 444.107 9.284 0.0010988768
249 2.396 NA NA NA -1.017 0.186 0.006 0.396 0.289 NA NA 7 -214.039 444.115 9.291 0.0010946917
457 2.416 NA NA NA -0.995 NA NA 0.361 0.312 -0.119 NA 6 -215.317 444.135 9.311 0.0010839898
649 2.279 NA NA NA -0.727 NA NA NA 0.282 NA 0.286 5 -216.557 444.166 9.342 0.0010672474
11 2.446 NA 0.242 NA -0.679 NA NA NA NA NA NA 4 -217.863 444.416 9.592 0.0009415836
203 2.417 NA 0.094 NA -0.948 NA NA 0.321 0.224 NA NA 6 -215.480 444.460 9.636 0.0009211244
202 2.356 0.058 NA NA -0.972 NA NA 0.323 0.266 NA NA 6 -215.538 444.576 9.753 0.0008691587
75 2.434 NA 0.216 NA -0.913 NA NA 0.304 NA NA NA 5 -216.763 444.578 9.754 0.0008685717
139 2.430 NA 0.131 NA -0.698 NA NA NA 0.213 NA NA 5 -216.764 444.582 9.758 0.0008668746
154 2.290 0.112 NA NA -0.750 0.169 NA NA 0.285 NA NA 6 -215.549 444.598 9.774 0.0008599298
889 1.896 NA NA NA -1.412 0.196 -0.159 0.839 NA 0.135 0.943 8 -212.972 444.610 9.787 0.0008544811

# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.chin.4 <- get.models(model.set.chin, subset = delta < 4)
avg_model.chin.4 <- model.avg(model.set.chin.4)
summary(avg_model.chin.4)
## 
## Call:
## model.avg(object = model.set.chin.4)
## 
## Component model call: 
## glm.nb(formula = abundance ~ <13 unique rhs>, data = b.chin.hab, 
##      na.action = na.fail, init.theta = <13 unique values>, link = log)
## 
## Component models: 
##              df  logLik   AICc delta weight
## 4/5/7/8/10    7 -209.39 434.82  0.00   0.21
## 4/7/8/10      6 -210.90 435.31  0.48   0.16
## 4/7/8/9/10    7 -210.39 436.81  1.99   0.08
## 4/5/6/7/8/10  8 -209.08 436.82  2.00   0.08
## 4/5/7/8/9/10  8 -209.26 437.18  2.36   0.06
## 3/4/5/7/8/10  8 -209.28 437.23  2.41   0.06
## 4/6/7/8/10    7 -210.67 437.39  2.56   0.06
## 2/4/7/8/10    7 -210.68 437.40  2.58   0.06
## 2/4/5/7/8/10  8 -209.38 437.42  2.60   0.06
## 1/4/5/7/8/10  8 -209.38 437.43  2.60   0.06
## 3/4/7/8/10    7 -210.76 437.56  2.74   0.05
## 1/4/7/8/10    7 -210.85 437.73  2.90   0.05
## 2/4/7/10      6 -212.64 438.79  3.96   0.03
## 
## Term codes: 
##       j2 meanturb     s.do s.J.date     s.pH    s.sal   s.temp stnwidth 
##        1        2        3        4        5        6        7        8 
##  vegelev     Year 
##        9       10 
## 
## Model-averaged coefficients:  
## (full average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)  1.833302   0.201193    0.205535   8.920  < 2e-16 ***
## s.J.date    -1.436520   0.249962    0.255335   5.626  < 2e-16 ***
## s.pH         0.090995   0.126463    0.127902   0.711  0.47681    
## s.temp       0.895418   0.270309    0.276154   3.242  0.00119 ** 
## stnwidth     0.310198   0.133734    0.136121   2.279  0.02268 *  
## Year         0.990602   0.303815    0.310332   3.192  0.00141 ** 
## vegelev     -0.014297   0.061036    0.061875   0.231  0.81727    
## s.sal       -0.015460   0.066537    0.067484   0.229  0.81880    
## s.do        -0.000654   0.048585    0.049450   0.013  0.98945    
## meanturb     0.013482   0.067085    0.067905   0.199  0.84262    
## j2          -0.002526   0.031215    0.031854   0.079  0.93680    
##  
## (conditional average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)  1.833302   0.201193    0.205535   8.920  < 2e-16 ***
## s.J.date    -1.436520   0.249962    0.255335   5.626  < 2e-16 ***
## s.pH         0.175321   0.126608    0.129364   1.355  0.17534    
## s.temp       0.895418   0.270309    0.276154   3.242  0.00119 ** 
## stnwidth     0.319258   0.124558    0.127191   2.510  0.01207 *  
## Year         0.990602   0.303815    0.310332   3.192  0.00141 ** 
## vegelev     -0.102398   0.132895    0.135646   0.755  0.45031    
## s.sal       -0.116291   0.146887    0.150102   0.775  0.43849    
## s.do        -0.005729   0.143703    0.146264   0.039  0.96876    
## meanturb     0.095407   0.155018    0.157524   0.606  0.54474    
## j2          -0.024218   0.093901    0.095938   0.252  0.80070    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      s.J.date s.temp Year stnwidth s.pH meanturb vegelev
## Importance:          1.00     1.00   1.00 0.97     0.52 0.14     0.14   
## N containing models:   13       13     13   12        6    3        2   
##                      s.sal s.do j2  
## Importance:          0.13  0.11 0.10
## N containing models:    2     2    2
chin.ci <- data.frame(confint(avg_model.chin.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4
chin.4.Rsq <- rsquared(model.set.chin.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.chin <- data.frame(avg_model.chin.4$msTable)
avg_model_components4.chin <- cbind(chin.4.Rsq, avg_model_4df.chin)
r = data.frame(Coeff = rownames(avg_model_4df.chin, rep(NA, length(avg_model_components4.chin))))
avg_model_components4.chin <- cbind(avg_model_components4.chin, r)
avg_model_components4.chin <- avg_model_components4.chin[, -c(6, 7)]
# write.csv(avg_model_components4.chin,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachchin.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

## Warning: package 'cowplot' was built under R version 3.4.3
## 
## Attaching package: 'cowplot'
## The following objects are masked from 'package:sjPlot':
## 
##     plot_grid, save_plot
## The following object is masked from 'package:ggplot2':
## 
##     ggsave

Check fit of top model by avg – including all parameters with weight >0.5

avg.chin <- glm.nb(abundance ~ s.J.date + Year + s.temp + stnwidth + s.pH, data = b.chin.hab, 
    na.action = "na.fail")

summary(avg.chin)  # AIC 432.9, deviance explained = 100*((129.963-64.196)/129.963) = 50.6 %
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + Year + s.temp + stnwidth + 
##     s.pH, data = b.chin.hab, na.action = "na.fail", init.theta = 1.414591464, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5158  -0.9501  -0.4285   0.3706   1.7588  
## 
## Coefficients:
##             Estimate Std. Error z value     Pr(>|z|)    
## (Intercept)   1.8395     0.1966   9.357      < 2e-16 ***
## s.J.date     -1.4488     0.2459  -5.892 0.0000000038 ***
## Year          0.9591     0.2963   3.237     0.001210 ** 
## s.temp        0.9078     0.2669   3.401     0.000672 ***
## stnwidth      0.3231     0.1185   2.725     0.006427 ** 
## s.pH          0.1744     0.1245   1.401     0.161140    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.4146) family taken to be 1)
## 
##     Null deviance: 129.963  on 62  degrees of freedom
## Residual deviance:  64.196  on 57  degrees of freedom
## AIC: 432.79
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.415 
##           Std. Err.:  0.270 
## 
##  2 x log-likelihood:  -418.787
vif(avg.chin)  # - Jday and temp still above 5 (5.055 and 5.93, respectively)
## s.J.date     Year   s.temp stnwidth     s.pH 
## 5.055297 1.639788 5.930079 1.067749 1.050073
### Plot residuals vs. fitted values
plot(fitted(avg.chin), resid(avg.chin), main = "Averaged Chinook GLM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = b.chin.hab$abundance, y = resid(avg.chin), main = "Q-Q Averaged Chinook GLM")
qqline(resid(avg.chin), col = 2)

Chum beach seine model

response = average abundance of chum salmon [per beach seine site]
This is 2nd of 4 beach models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from beach seine data

### Separate into groups; automatically removes 3 x unidentified larval and 28
### x '0' abundance rows 2
b.2 <- b[which(b$Species == "Chum"), ]

### standardize continuous variables
b.2$s.temp <- standardize(b.2$Temp.surf, centerFun = mean, scaleFun = sd)
b.2$s.sal <- standardize(b.2$Sal.surf, centerFun = mean, scaleFun = sd)
b.2$s.do <- standardize(b.2$DOmg.surf, centerFun = mean, scaleFun = sd)
b.2$s.pH <- standardize(b.2$pH.surf, centerFun = mean, scaleFun = sd)
b.2$s.J.date <- standardize(b.2$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.2$j2 <- b.2$s.J.date^2

# 1 summarize by site-day
b.chum <- ddply(b.2, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.chum)

Habitat variables

## Add marsh habitat variables for each site
b.chum.hab <- b.chum
b.chum.hab$stnwidth <- hab2[match(b.chum.hab$Site, hab$Site), 2]
b.chum.hab$vegelev <- hab2[match(b.chum.hab$Site, hab$Site), 3]
b.chum.hab$shtdensity <- hab2[match(b.chum.hab$Site, hab$Site), 4]
b.chum.hab$shtdenhigh <- hab2[match(b.chum.hab$Site, hab$Site), 5]
b.chum.hab$tcelev <- hab2[match(b.chum.hab$Site, hab$Site), 6]
b.chum.hab$angbank <- hab2[match(b.chum.hab$Site, hab$Site), 7]
b.chum.hab$meanturb <- hab2[match(b.chum.hab$Site, hab$Site), 8]

## standardize habitat variables using 'standardize' function from robustHD
## package
b.chum.hab$stnwidth <- standardize(b.chum.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.chum.hab$stnwidth <- as.numeric(b.chum.hab$stnwidth)
b.chum.hab$vegelev <- standardize(b.chum.hab$vegelev, centerFun = mean, scaleFun = sd)
b.chum.hab$vegelev <- as.numeric(b.chum.hab$vegelev)
b.chum.hab$shtdensity <- standardize(b.chum.hab$shtdensity, centerFun = mean, 
    scaleFun = sd)
b.chum.hab$shtdensity <- as.numeric(b.chum.hab$shtdensity)
b.chum.hab$shtdenhigh <- standardize(b.chum.hab$shtdenhigh, centerFun = mean, 
    scaleFun = sd)
b.chum.hab$shtdenhigh <- as.numeric(b.chum.hab$shtdenhigh)
b.chum.hab$tcelev <- standardize(b.chum.hab$tcelev, centerFun = mean, scaleFun = sd)
b.chum.hab$tcelev <- as.numeric(b.chum.hab$tcelev)
b.chum.hab$angbank <- standardize(b.chum.hab$angbank, centerFun = mean, scaleFun = sd)
b.chum.hab$angbank <- as.numeric(b.chum.hab$angbank)
b.chum.hab$meanturb <- standardize(b.chum.hab$meanturb, centerFun = mean, scaleFun = sd)
b.chum.hab$meanturb <- as.numeric(b.chum.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors


# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.chum.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
##     stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank + 
##     meanturb
## 
## Complete :
##          (Intercept)          s.J.date             j2                  
## tcelev                      0                    0                    0
## angbank                     0                    0                    0
## meanturb                    0                    0                    0
##          Year                 s.temp               s.sal               
## tcelev                      0                    0                    0
## angbank                     0                    0                    0
## meanturb                    0                    0                    0
##          s.do                 s.pH                 stnwidth            
## tcelev                      0                    0           8871/14840
## angbank                     0                    0       -133424/254871
## meanturb                    0                    0          27432/68227
##          vegelev              shtdenhigh           shtdensity          
## tcelev   605155061/1293820358       -376609/756041            1235/3802
## angbank             2957/4166          -9279/22789          -3094/47891
## meanturb     -3080706/3772433      -871961/1287673         -35732/62255

# We determined hab variables for Chum should be reduced to meanturb, and
# vegelev based on significance and collinearity when doing model
# exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + vegelev + 
    meanturb, data = b.chum.hab))  ##all less than 5
## s.J.date       j2     Year   s.temp    s.sal     s.do     s.pH  vegelev 
## 2.603210 2.909987 2.583522 3.466341 3.020429 3.369835 1.630916 2.041842 
## meanturb 
## 1.165110


## Pearson Corr with all vars
year <- b.chum.hab$Year
jday <- b.chum.hab$s.J.date
j2 <- b.chum.hab$j2
temp <- b.chum.hab$s.temp
do <- b.chum.hab$s.do
ph <- b.chum.hab$s.pH
sal <- b.chum.hab$s.sal
veg <- b.chum.hab$vegelev
turb <- b.chum.hab$meanturb
width <- b.chum.hab$stnwidth
shootden <- b.chum.hab$shtdensity
shtdenhigh <- b.chum.hab$shtdenhigh
elev <- b.chum.hab$tcelev
angbank <- b.chum.hab$angbank

habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width, 
    shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine chum salmon dataset")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chummarsh.pdf")

Model selection

Full model

## Chum: full model with top habitat variables
chum <- glm.nb(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
    vegelev + meanturb, data = b.chum.hab, na.action = "na.fail")
summary(chum)  #AIC 233.64
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.temp + 
##     s.sal + s.do + s.pH + vegelev + meanturb, data = b.chum.hab, 
##     na.action = "na.fail", init.theta = 0.9805681787, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8251  -1.2597  -0.3236   0.5022   1.4294  
## 
## Coefficients:
##             Estimate Std. Error z value    Pr(>|z|)    
## (Intercept)   2.2029     0.4342   5.073 0.000000391 ***
## s.J.date     -0.5285     0.2968  -1.781      0.0749 .  
## j2           -0.1948     0.2205  -0.884      0.3770    
## Year          0.9408     0.6445   1.460      0.1443    
## s.temp        0.1084     0.3458   0.314      0.7539    
## s.sal        -0.4729     0.3153  -1.500      0.1336    
## s.do          0.5665     0.3698   1.532      0.1256    
## s.pH          0.3145     0.2707   1.162      0.2452    
## vegelev       0.3977     0.2829   1.406      0.1598    
## meanturb     -0.3448     0.2157  -1.599      0.1098    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9806) family taken to be 1)
## 
##     Null deviance: 58.891  on 30  degrees of freedom
## Residual deviance: 33.073  on 21  degrees of freedom
## AIC: 233.64
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.981 
##           Std. Err.:  0.255 
## 
##  2 x log-likelihood:  -211.642
## Now look at the variance inflation factors of this full model
vif(chum)
## s.J.date       j2     Year   s.temp    s.sal     s.do     s.pH  vegelev 
## 2.366396 2.638545 2.607406 3.100914 3.139654 3.062227 1.634527 2.029153 
## meanturb 
## 1.210263
# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(chum, type = "ma")
## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.
## No outliers detected.
## NULL

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 0.12753
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.92707
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 2.3502e-16
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 0.12753
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 0.92707
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 2.3502e-16
## --------------------
## Check significance of terms when they entered the model...
## Anova:
## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'

## Analysis of Deviance Table
## 
## Model: Negative Binomial(0.9806), link: log
## 
## Response: abundance
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        30     58.891              
## s.J.date  1   0.9566        29     57.935 0.3280324    
## j2        1  14.0090        28     43.926 0.0001819 ***
## Year      1   3.6769        27     40.249 0.0551719 .  
## s.temp    1   0.0207        26     40.228 0.8855816    
## s.sal     1   0.5239        25     39.704 0.4691742    
## s.do      1   2.8477        24     36.856 0.0915062 .  
## s.pH      1   0.2736        23     36.583 0.6009037    
## vegelev   1   1.0859        22     35.497 0.2973723    
## meanturb  1   2.4243        21     33.073 0.1194640    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# See pseudo R2 for full model
rsquared(chum, aicc = TRUE)
##    Class            Family Link  n R.squared
## 1 negbin Negative Binomial  log 31 0.4384135

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (chum has 31 observations). Note that due to the small sample size, even the top model has a very low weight (0.08), and consequently many models (42) are included in the averaging set due to very small changes in delta AICc with different iterations.

# Generate model set
model.set.chum <- dredge(chum)

# Create model selection table
model_table.chum <- model.sel(model.set.chum)
options(scipen = 7)
names(model_table.chum) <- c("(Int)", "Jday^2", "Mean turb", "DO", "Jday", "pH", 
    "Sal", "Temp", "Veg. elev", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(model_table.chum, n = 100), digits = 3)
(Int) Jday^2 Mean turb DO Jday pH Sal Temp Veg. elev Year df LL AICc delta weight
13 2.621 NA NA 0.923 -0.367 NA NA NA NA NA 4 -108.612 226.763 0.000 0.046856556
5 2.646 NA NA 0.838 NA NA NA NA NA NA 3 -110.013 226.914 0.151 0.043447042
69 2.599 NA NA 0.886 NA NA NA -0.319 NA NA 4 -108.998 227.534 0.771 0.031868966
15 2.572 NA -0.302 0.917 -0.377 NA NA NA NA NA 5 -107.587 227.574 0.811 0.031233524
7 2.610 NA -0.290 0.847 NA NA NA NA NA NA 4 -109.159 227.857 1.094 0.027115205
269 2.306 NA NA 0.900 -0.411 NA NA NA NA 0.474 5 -108.047 228.495 1.732 0.019709864
133 2.635 NA NA 0.885 NA NA NA NA 0.201 NA 4 -109.623 228.784 2.021 0.017059660
37 2.610 NA NA 0.831 NA NA 0.189 NA NA NA 4 -109.686 228.911 2.149 0.016003777
77 2.592 NA NA 0.920 -0.280 NA NA -0.207 NA NA 5 -108.298 228.996 2.233 0.015341997
261 2.433 NA NA 0.826 NA NA NA NA NA 0.329 4 -109.739 229.017 2.254 0.015183630
6 2.806 -0.173 NA 0.629 NA NA NA NA NA NA 4 -109.785 229.108 2.345 0.014503333
2 3.170 -0.516 NA NA NA NA NA NA NA NA 3 -111.114 229.117 2.354 0.014440692
71 2.575 NA -0.240 0.881 NA NA NA -0.269 NA NA 5 -108.423 229.245 2.482 0.013543770
141 2.611 NA NA 0.944 -0.340 NA NA NA 0.139 NA 5 -108.423 229.246 2.484 0.013535084
18 3.109 -0.473 NA NA NA 0.313 NA NA NA NA 4 -109.861 229.261 2.498 0.013439000
66 3.116 -0.540 NA NA NA NA NA -0.376 NA NA 4 -109.873 229.283 2.521 0.013286939
14 2.735 -0.128 NA 0.767 -0.352 NA NA NA NA NA 5 -108.466 229.332 2.569 0.012967934
21 2.649 NA NA 0.796 NA 0.075 NA NA NA NA 4 -109.958 229.455 2.692 0.012196116
45 2.603 NA NA 0.908 -0.342 NA 0.091 NA NA NA 5 -108.540 229.480 2.717 0.012043351
70 2.796 -0.225 NA 0.612 NA NA NA -0.345 NA NA 5 -108.576 229.551 2.788 0.011623464
29 2.623 NA NA 0.894 -0.365 0.046 NA NA NA NA 5 -108.588 229.576 2.814 0.011476432
135 2.595 NA -0.303 0.901 NA NA NA NA 0.221 NA 5 -108.647 229.694 2.931 0.010821136
197 2.589 NA NA 0.923 NA NA NA -0.308 0.173 NA 5 -108.672 229.743 2.980 0.010558866
271 2.301 NA -0.289 0.896 -0.417 NA NA NA NA 0.416 6 -107.127 229.754 2.991 0.010504088
258 2.733 -0.522 NA NA NA NA NA NA NA 0.654 4 -110.114 229.767 3.004 0.010435365
20 3.045 -0.446 -0.340 NA NA 0.440 NA NA NA NA 5 -108.724 229.847 3.084 0.010024106
82 3.070 -0.501 NA NA NA 0.292 NA -0.336 NA NA 5 -108.759 229.917 3.154 0.009678037
39 2.572 NA -0.292 0.842 NA NA 0.193 NA NA NA 5 -108.785 229.971 3.208 0.009423040
143 2.559 NA -0.306 0.945 -0.347 NA NA NA 0.153 NA 6 -107.340 230.180 3.417 0.008486253
23 2.611 NA -0.332 0.746 NA 0.186 NA NA NA NA 5 -108.898 230.196 3.433 0.008419450
31 2.574 NA -0.340 0.821 -0.374 0.157 NA NA NA NA 6 -107.367 230.234 3.471 0.008260439
325 2.685 NA NA 0.898 NA NA NA -0.368 NA -0.144 5 -108.967 230.334 3.572 0.007856236
263 2.440 NA -0.279 0.838 NA NA NA NA NA 0.264 5 -108.975 230.350 3.587 0.007796556
85 2.601 NA NA 0.868 NA 0.031 NA -0.315 NA NA 5 -108.989 230.378 3.615 0.007687679
101 2.596 NA NA 0.883 NA NA 0.027 -0.305 NA NA 5 -108.993 230.385 3.623 0.007658219
79 2.558 NA -0.281 0.914 -0.323 NA NA -0.133 NA NA 6 -107.453 230.405 3.642 0.007582632
8 2.734 -0.133 -0.276 0.688 NA NA NA NA NA NA 5 -109.024 230.447 3.684 0.007425643
47 2.553 NA -0.302 0.904 -0.352 NA 0.095 NA NA NA 6 -107.499 230.498 3.735 0.007239276
266 2.603 -0.523 NA NA -0.326 NA NA NA NA 0.802 5 -109.087 230.573 3.811 0.006971391
26 3.097 -0.485 NA NA -0.269 0.316 NA NA NA NA 5 -109.096 230.591 3.828 0.006909639
16 2.634 -0.069 -0.291 0.836 -0.367 NA NA NA NA NA 6 -107.546 230.593 3.830 0.006903870
10 3.161 -0.526 NA NA -0.242 NA NA NA NA NA 4 -110.546 230.631 3.868 0.006773211
4 3.142 -0.512 -0.219 NA NA NA NA NA NA NA 4 -110.648 230.834 4.071 0.006119945
34 3.101 -0.498 NA NA NA NA 0.236 NA NA NA 4 -110.711 230.960 4.197 0.005745496
270 2.423 -0.194 NA 0.658 -0.389 NA NA NA NA 0.554 6 -107.731 230.961 4.199 0.005742065
262 2.581 -0.244 NA 0.522 NA NA NA NA NA 0.446 5 -109.324 231.047 4.284 0.005500541
274 2.808 -0.484 NA NA NA 0.265 NA NA NA 0.467 5 -109.340 231.081 4.318 0.005409451
389 2.454 NA NA 0.870 NA NA NA NA 0.178 0.281 5 -109.424 231.249 4.486 0.004973398
28 3.015 -0.450 -0.342 NA -0.280 0.438 NA NA NA NA 6 -107.875 231.251 4.488 0.004968677
38 2.770 -0.177 NA 0.615 NA NA 0.195 NA NA NA 5 -109.438 231.277 4.514 0.004904452
149 2.637 NA NA 0.817 NA 0.124 NA NA 0.230 NA 5 -109.473 231.347 4.584 0.004735403
397 2.318 NA NA 0.919 -0.389 NA NA NA 0.105 0.445 6 -107.934 231.367 4.605 0.004686947
22 2.903 -0.272 NA 0.404 NA 0.186 NA NA NA NA 5 -109.514 231.428 4.666 0.004546139
134 2.736 -0.108 NA 0.747 NA NA NA NA 0.169 NA 5 -109.539 231.478 4.715 0.004435294
301 2.247 NA NA 0.909 -0.447 NA -0.099 NA NA 0.591 6 -107.992 231.483 4.720 0.004423265
199 2.563 NA -0.249 0.922 NA NA NA -0.248 0.187 NA 6 -108.027 231.554 4.791 0.004270692
165 2.620 NA NA 0.867 NA NA 0.093 NA 0.140 NA 5 -109.579 231.558 4.795 0.004261933
333 2.278 NA NA 0.899 -0.431 NA NA 0.038 NA 0.524 6 -108.043 231.587 4.824 0.004199797
285 2.309 NA NA 0.895 -0.410 0.008 NA NA NA 0.471 6 -108.047 231.594 4.831 0.004185970
78 2.732 -0.163 NA 0.721 -0.243 NA NA -0.234 NA NA 6 -108.071 231.643 4.880 0.004084400
293 2.502 NA NA 0.828 NA NA 0.134 NA NA 0.183 5 -109.626 231.652 4.889 0.004065417
68 3.098 -0.533 -0.150 NA NA NA NA -0.340 NA NA 5 -109.658 231.716 4.953 0.003937852
322 2.894 -0.538 NA NA NA NA NA -0.265 NA 0.361 5 -109.663 231.726 4.963 0.003917347
84 3.026 -0.472 -0.261 NA NA 0.382 NA -0.255 NA NA 6 -108.117 231.734 4.971 0.003902126
50 3.071 -0.465 NA NA NA 0.292 0.152 NA NA NA 5 -109.669 231.737 4.974 0.003895890
53 2.613 NA NA 0.807 NA 0.043 0.180 NA NA NA 5 -109.669 231.738 4.975 0.003894211
205 2.583 NA NA 0.942 -0.253 NA NA -0.204 0.132 NA 6 -108.120 231.739 4.976 0.003891784
130 3.170 -0.516 NA NA NA NA NA NA -0.022 NA 4 -111.109 231.757 4.994 0.003856916
72 2.740 -0.184 -0.210 0.659 NA NA NA -0.294 NA NA 6 -108.148 231.796 5.033 0.003783293
151 2.591 NA -0.353 0.767 NA 0.246 NA NA 0.265 NA 6 -108.170 231.840 5.077 0.003700746
277 2.446 NA NA 0.802 NA 0.043 NA NA NA 0.310 5 -109.723 231.846 5.083 0.003690078
146 3.094 -0.462 NA NA NA 0.345 NA NA 0.127 NA 5 -109.724 231.848 5.085 0.003686308
260 2.735 -0.520 -0.199 NA NA NA NA NA NA 0.616 5 -109.736 231.873 5.110 0.003640283
282 2.677 -0.489 NA NA -0.336 0.263 NA NA NA 0.620 6 -108.192 231.883 5.120 0.003621802
74 3.111 -0.537 NA NA -0.091 NA NA -0.337 NA NA 5 -109.806 232.013 5.250 0.003394477
93 2.594 NA NA 0.898 -0.281 0.035 NA -0.203 NA NA 6 -108.285 232.070 5.307 0.003299132
194 3.120 -0.543 NA NA NA NA NA -0.378 -0.049 NA 5 -109.847 232.095 5.332 0.003257983
109 2.592 NA NA 0.919 -0.280 NA 0.004 -0.205 NA NA 6 -108.298 232.096 5.333 0.003256704
30 2.803 -0.199 NA 0.600 -0.337 0.131 NA NA NA NA 6 -108.316 232.133 5.370 0.003196857
87 2.579 NA -0.268 0.811 NA 0.118 NA -0.245 NA NA 6 -108.318 232.137 5.374 0.003190311
98 3.112 -0.538 NA NA NA NA 0.018 -0.367 NA NA 5 -109.871 232.142 5.379 0.003182667
157 2.613 NA NA 0.889 -0.333 0.093 NA NA 0.170 NA 6 -108.332 232.165 5.402 0.003146370
142 2.691 -0.089 NA 0.832 -0.335 NA NA NA 0.111 NA 6 -108.360 232.219 5.456 0.003061326
86 2.873 -0.303 NA 0.420 NA 0.159 NA -0.333 NA NA 6 -108.380 232.260 5.497 0.003000189
327 2.666 NA -0.239 0.892 NA NA NA -0.321 NA -0.152 6 -108.388 232.277 5.514 0.002975030
46 2.717 -0.129 NA 0.750 -0.326 NA 0.093 NA NA NA 6 -108.391 232.283 5.520 0.002966136
103 2.567 NA -0.246 0.875 NA NA 0.065 -0.233 NA NA 6 -108.393 232.286 5.523 0.002961601
173 2.612 NA NA 0.946 -0.341 NA -0.004 NA 0.142 NA 6 -108.423 232.346 5.583 0.002873089
148 3.017 -0.429 -0.356 NA NA 0.489 NA NA 0.172 NA 6 -108.437 232.373 5.610 0.002834510
198 2.744 -0.176 NA 0.697 NA NA NA -0.332 0.119 NA 6 -108.438 232.377 5.614 0.002829707
229 2.612 NA NA 0.988 NA NA -0.259 -0.425 0.318 NA 6 -108.439 232.378 5.615 0.002828478
24 2.864 -0.269 -0.330 0.359 NA 0.307 NA NA NA NA 6 -108.441 232.381 5.619 0.002822976
276 2.845 -0.458 -0.304 NA NA 0.383 NA NA NA 0.323 6 -108.473 232.446 5.683 0.002733817
61 2.606 NA NA 0.885 -0.343 0.038 0.086 NA NA NA 6 -108.524 232.548 5.785 0.002597375
391 2.461 NA -0.294 0.891 NA NA NA NA 0.205 0.211 6 -108.529 232.558 5.795 0.002584493
12 3.123 -0.518 -0.214 NA -0.242 NA NA NA NA NA 5 -110.083 232.567 5.804 0.002573197
386 2.733 -0.523 NA NA NA NA NA NA -0.031 0.655 5 -110.104 232.608 5.845 0.002520464
90 3.060 -0.495 NA NA -0.152 0.301 NA -0.272 NA NA 6 -108.557 232.615 5.852 0.002512494
290 2.745 -0.520 NA NA NA NA 0.026 NA NA 0.625 5 -110.111 232.622 5.859 0.002503511
102 2.795 -0.225 NA 0.612 NA NA 0.006 -0.342 NA NA 6 -108.575 232.651 5.888 0.002467574

# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.chum.4 <- get.models(model.set.chum, subset = delta < 4)
avg_model.chum.4 <- model.avg(model.set.chum.4)
summary(avg_model.chum.4)
## 
## Call:
## model.avg(object = model.set.chum.4)
## 
## Component model call: 
## glm.nb(formula = abundance ~ <42 unique rhs>, data = b.chum.hab, 
##      na.action = na.fail, init.theta = <42 unique values>, link = log)
## 
## Component models: 
##      df  logLik   AICc delta weight
## 34    4 -108.61 226.76  0.00   0.08
## 3     3 -110.01 226.91  0.15   0.07
## 37    4 -109.00 227.53  0.77   0.05
## 234   5 -107.59 227.57  0.81   0.05
## 23    4 -109.16 227.86  1.09   0.05
## 349   5 -108.05 228.49  1.73   0.03
## 38    4 -109.62 228.78  2.02   0.03
## 36    4 -109.69 228.91  2.15   0.03
## 347   5 -108.30 229.00  2.23   0.03
## 39    4 -109.74 229.02  2.25   0.03
## 13    4 -109.78 229.11  2.35   0.02
## 1     3 -111.11 229.12  2.35   0.02
## 237   5 -108.42 229.25  2.48   0.02
## 348   5 -108.42 229.25  2.48   0.02
## 15    4 -109.86 229.26  2.50   0.02
## 17    4 -109.87 229.28  2.52   0.02
## 134   5 -108.47 229.33  2.57   0.02
## 35    4 -109.96 229.45  2.69   0.02
## 346   5 -108.54 229.48  2.72   0.02
## 137   5 -108.58 229.55  2.79   0.02
## 345   5 -108.59 229.58  2.81   0.02
## 238   5 -108.65 229.69  2.93   0.02
## 378   5 -108.67 229.74  2.98   0.02
## 2349  6 -107.13 229.75  2.99   0.02
## 19    4 -110.11 229.77  3.00   0.02
## 125   5 -108.72 229.85  3.08   0.02
## 157   5 -108.76 229.92  3.15   0.02
## 236   5 -108.79 229.97  3.21   0.02
## 2348  6 -107.34 230.18  3.42   0.01
## 235   5 -108.90 230.20  3.43   0.01
## 2345  6 -107.37 230.23  3.47   0.01
## 379   5 -108.97 230.33  3.57   0.01
## 239   5 -108.97 230.35  3.59   0.01
## 357   5 -108.99 230.38  3.61   0.01
## 367   5 -108.99 230.39  3.62   0.01
## 2347  6 -107.45 230.41  3.64   0.01
## 123   5 -109.02 230.45  3.68   0.01
## 2346  6 -107.50 230.50  3.74   0.01
## 149   5 -109.09 230.57  3.81   0.01
## 145   5 -109.10 230.59  3.83   0.01
## 1234  6 -107.55 230.59  3.83   0.01
## 14    4 -110.55 230.63  3.87   0.01
## 
## Term codes: 
##       j2 meanturb     s.do s.J.date     s.pH    s.sal   s.temp  vegelev 
##        1        2        3        4        5        6        7        8 
##     Year 
##        9 
## 
## Model-averaged coefficients:  
## (full average) 
##             Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)  2.66487    0.31503     0.32401   8.225   <2e-16 ***
## s.do         0.72465    0.39016     0.39608   1.830   0.0673 .  
## s.J.date    -0.14133    0.22008     0.22390   0.631   0.5279    
## s.temp      -0.07041    0.17154     0.17492   0.403   0.6873    
## meanturb    -0.08859    0.17905     0.18265   0.485   0.6277    
## Year         0.05483    0.23053     0.23587   0.232   0.8162    
## vegelev      0.01855    0.08886     0.09141   0.203   0.8392    
## s.sal        0.01168    0.07314     0.07546   0.155   0.8770    
## j2          -0.09338    0.20253     0.20422   0.457   0.6475    
## s.pH         0.03107    0.13235     0.13566   0.229   0.8188    
##  
## (conditional average) 
##             Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   2.6649     0.3150      0.3240   8.225  < 2e-16 ***
## s.do          0.8600     0.2535      0.2642   3.255  0.00113 ** 
## s.J.date     -0.3547     0.2142      0.2239   1.584  0.11314    
## s.temp       -0.3010     0.2375      0.2478   1.215  0.22445    
## meanturb     -0.2961     0.2138      0.2237   1.324  0.18561    
## Year          0.4084     0.5014      0.5196   0.786  0.43190    
## vegelev       0.1792     0.2179      0.2279   0.786  0.43162    
## s.sal         0.1303     0.2103      0.2193   0.594  0.55234    
## j2           -0.3753     0.2430      0.2486   1.509  0.13117    
## s.pH          0.2061     0.2831      0.2933   0.703  0.48223    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      s.do s.J.date meanturb j2   s.temp s.pH Year vegelev
## Importance:          0.84 0.40     0.30     0.25 0.23   0.15 0.13 0.10   
## N containing models:   33   17       15       14   11      9    7    5   
##                      s.sal
## Importance:          0.09 
## N containing models:    5
chum.ci <- data.frame(confint(avg_model.chum.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4
chum.4.Rsq <- rsquared(model.set.chum.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.chum <- data.frame(avg_model.chum.4$msTable)
avg_model_components4.chum <- cbind(chum.4.Rsq, avg_model_4df.chum)
r = data.frame(Coeff = rownames(avg_model_4df.chum, rep(NA, length(avg_model_components4.chum))))
avg_model_components4.chum <- cbind(avg_model_components4.chum, r)
avg_model_components4.chum <- avg_model_components4.chum[, -c(6, 7)]
# write.csv(avg_model_components4.chum,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachchum.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

Check fit of top model by avg – including all parameters with weight >0.5

avg.chum <- glm.nb(abundance ~ s.do, data = b.chum.hab, na.action = "na.fail")

summary(avg.chum)  # AIC 226.03  
## 
## Call:
## glm.nb(formula = abundance ~ s.do, data = b.chum.hab, na.action = "na.fail", 
##     init.theta = 0.7611315465, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7338  -1.2327  -0.4714   0.4236   1.7209  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.6465     0.2152  12.298  < 2e-16 ***
## s.do          0.8380     0.2296   3.651 0.000262 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.7611) family taken to be 1)
## 
##     Null deviance: 46.900  on 30  degrees of freedom
## Residual deviance: 34.031  on 29  degrees of freedom
## AIC: 226.03
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.761 
##           Std. Err.:  0.186 
## 
##  2 x log-likelihood:  -220.025
# vif(avg.chum) ## Can't run VIF because only one variable!

### Plot residuals vs. fitted values
plot(fitted(avg.chum), resid(avg.chum), main = "Averaged Chum GLM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = b.chum.hab$abundance, y = resid(avg.chum), main = "Q-Q Averaged Chum GLM")
qqline(resid(avg.chum), col = 2)

Migratory beach seine model

response = average abundance of migratory fishes [per beach seine site]
This is 3rd of 4 beach models for (1)Chinook, (2)migra, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from beach seine data

### Grab just migratory catch - automatically removes any 0s or unidentified
### species.  3
b.3 <- b[which(b$class == "migratory"), ]
b.3 <- b.3[which(b.3$Species != "Chinook"), ]
b.3 <- b.3[which(b.3$Species != "Chum"), ]

### standardize continuous variables
b.3$s.temp <- standardize(b.3$Temp.surf, centerFun = mean, scaleFun = sd)
b.3$s.sal <- standardize(b.3$Sal.surf, centerFun = mean, scaleFun = sd)
b.3$s.do <- standardize(b.3$DOmg.surf, centerFun = mean, scaleFun = sd)
b.3$s.pH <- standardize(b.3$pH.surf, centerFun = mean, scaleFun = sd)
b.3$s.J.date <- standardize(b.3$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.3$j2 <- b.3$s.J.date^2

# 1 summarize by site-day
b.migra <- ddply(b.3, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.migra)

Habitat variables

## Add marsh habitat variables for each site
b.migra.hab <- b.migra
b.migra.hab$stnwidth <- hab2[match(b.migra.hab$Site, hab$Site), 2]
b.migra.hab$vegelev <- hab2[match(b.migra.hab$Site, hab$Site), 3]
b.migra.hab$shtdensity <- hab2[match(b.migra.hab$Site, hab$Site), 4]
b.migra.hab$shtdenhigh <- hab2[match(b.migra.hab$Site, hab$Site), 5]
b.migra.hab$tcelev <- hab2[match(b.migra.hab$Site, hab$Site), 6]
b.migra.hab$angbank <- hab2[match(b.migra.hab$Site, hab$Site), 7]
b.migra.hab$meanturb <- hab2[match(b.migra.hab$Site, hab$Site), 8]

## standardize habitat variables using 'standardize' function from robustHD
## package
b.migra.hab$stnwidth <- standardize(b.migra.hab$stnwidth, centerFun = mean, 
    scaleFun = sd)
b.migra.hab$stnwidth <- as.numeric(b.migra.hab$stnwidth)
b.migra.hab$vegelev <- standardize(b.migra.hab$vegelev, centerFun = mean, scaleFun = sd)
b.migra.hab$vegelev <- as.numeric(b.migra.hab$vegelev)
b.migra.hab$shtdensity <- standardize(b.migra.hab$shtdensity, centerFun = mean, 
    scaleFun = sd)
b.migra.hab$shtdensity <- as.numeric(b.migra.hab$shtdensity)
b.migra.hab$shtdenhigh <- standardize(b.migra.hab$shtdenhigh, centerFun = mean, 
    scaleFun = sd)
b.migra.hab$shtdenhigh <- as.numeric(b.migra.hab$shtdenhigh)
b.migra.hab$tcelev <- standardize(b.migra.hab$tcelev, centerFun = mean, scaleFun = sd)
b.migra.hab$tcelev <- as.numeric(b.migra.hab$tcelev)
b.migra.hab$angbank <- standardize(b.migra.hab$angbank, centerFun = mean, scaleFun = sd)
b.migra.hab$angbank <- as.numeric(b.migra.hab$angbank)
b.migra.hab$meanturb <- standardize(b.migra.hab$meanturb, centerFun = mean, 
    scaleFun = sd)
b.migra.hab$meanturb <- as.numeric(b.migra.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.migra.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
##     stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank + 
##     meanturb
## 
## Complete :
##          (Intercept)        s.J.date           j2                
## tcelev                    0                  0                  0
## angbank                   0                  0                  0
## meanturb                  0                  0                  0
##          Year               s.temp             s.sal             
## tcelev                    0                  0                  0
## angbank                   0                  0                  0
## meanturb                  0                  0                  0
##          s.do               s.pH               stnwidth          
## tcelev                    0                  0        12768/21383
## angbank                   0                  0         -1190/2477
## meanturb                  0                  0        35115/75149
##          vegelev            shtdenhigh         shtdensity        
## tcelev          26253/53636     -207653/439893    1598375/5117629
## angbank           6410/9393        -4595/12962      -15299/268071
## meanturb -12294629/12365285      -86771/116221       -43634/67957

# We determined hab variables for Migratory should be reduced to meanturb,
# and vegelev based on significance and collinearity when doing model
# exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + meanturb, data = b.migra.hab))  ##all less than 5
##  s.J.date        j2      Year    s.temp     s.sal      s.do      s.pH 
## 16.198637  6.067386 19.980166  6.905689 12.298055  7.072571  5.577524 
##  stnwidth   vegelev  meanturb 
##  2.771200  1.550012  3.061611


## Pearson Corr with all vars
year <- b.migra.hab$Year
jday <- b.migra.hab$s.J.date
j2 <- b.migra.hab$j2
temp <- b.migra.hab$s.temp
do <- b.migra.hab$s.do
ph <- b.migra.hab$s.pH
sal <- b.migra.hab$s.sal
veg <- b.migra.hab$vegelev
turb <- b.migra.hab$meanturb
width <- b.migra.hab$stnwidth
shootden <- b.migra.hab$shtdensity
shtdenhigh <- b.migra.hab$shtdenhigh
elev <- b.migra.hab$tcelev
angbank <- b.migra.hab$angbank

habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width, 
    shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine migratory fish dataset")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_migratorymarsh.pdf")

Model selection

Full model

Due to the very small sample size for other migratory fishes (# observations = 17), there were issues with the model assumptions with several parameters, and the global model was reduced to the maximum number of parameters than met assumptions and allowed the model to converge. Note that models with fewer, different parameters also converged (namely with temperature but without mean turb and stnwidth), but the AIC was worse and we felt the global model should include as many parameters as possible. The outliers in this data set caused a few parameters to always give errors: Jdate, J2 and pH. These were removed from the global model. Year had VIF over 10, so removed and AIC improved by delta 2 and VIF all were reduced to below 1.

## Migratory: full model with top habitat variables
migra <- glm.nb(abundance ~ s.do + s.sal + stnwidth + vegelev + meanturb, data = b.migra.hab, 
    na.action = "na.fail")
summary(migra)  #AIC 104.31
## 
## Call:
## glm.nb(formula = abundance ~ s.do + s.sal + stnwidth + vegelev + 
##     meanturb, data = b.migra.hab, na.action = "na.fail", init.theta = 1.175221356, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.1737  -0.9135  -0.2309   0.2802   1.4894  
## 
## Coefficients:
##             Estimate Std. Error z value      Pr(>|z|)    
## (Intercept)   1.6004     0.2689   5.952 0.00000000265 ***
## s.do          0.4460     0.3076   1.450       0.14708    
## s.sal        -0.8763     0.2669  -3.283       0.00103 ** 
## stnwidth      0.1927     0.3750   0.514       0.60726    
## vegelev       0.8562     0.2958   2.895       0.00380 ** 
## meanturb      0.4815     0.3838   1.255       0.20965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.1752) family taken to be 1)
## 
##     Null deviance: 44.223  on 16  degrees of freedom
## Residual deviance: 15.147  on 11  degrees of freedom
## AIC: 104.31
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.175 
##           Std. Err.:  0.431 
## 
##  2 x log-likelihood:  -90.310
vif(migra)
##     s.do    s.sal stnwidth  vegelev meanturb 
## 1.353410 1.194678 1.845305 1.218116 1.743007
# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glm(migra, type = "ma")
## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.
## 4 outliers removed in updated model.
## # A tibble: 2 x 2
##   models     aic
##   <chr>    <dbl>
## 1 original 104. 
## 2 updated   81.6

## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 1.0818
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 1.1311
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 3.9819e-17
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 1.0969
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used
## at 1.0818
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius
## 1.1311
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal
## condition number 3.9819e-17
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other
## near singularities as well. 1.0969
## --------------------
## Check significance of terms when they entered the model...
## Anova:
## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'

## Analysis of Deviance Table
## 
## Model: Negative Binomial(1.1752), link: log
## 
## Response: abundance
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev   Pr(>Chi)    
## NULL                        16     44.223               
## s.do      1   0.3035        15     43.920   0.581700    
## s.sal     1  16.3559        14     27.564 0.00005249 ***
## stnwidth  1   2.6652        13     24.899   0.102563    
## vegelev   1   8.3786        12     16.520   0.003797 ** 
## meanturb  1   1.3734        11     15.147   0.241235    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# See pseudo R2 for full model
rsquared(migra, aicc = TRUE)
##    Class            Family Link  n R.squared
## 1 negbin Negative Binomial  log 17  0.657492

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (migra has 17 observations).

# Generate model set
model.set.migra <- dredge(migra)

# Create model selection table
model_table.migra <- model.sel(model.set.migra)
options(scipen = 7)
names(model_table.migra) <- c("(Int)", "Mean turb", "DO", "Sal", "Width", "Veg. elev", 
    "df", "LL", "AICc", "delta", "weight")
kable(head(model_table.migra, n = 100), digits = 3)
(Int) Mean turb DO Sal Width Veg. elev df LL AICc delta weight
21 1.776 NA NA -0.791 NA 0.752 4 -47.440 106.212 0.000 0.29547923989
22 1.656 0.528 NA -0.757 NA 0.831 5 -46.113 107.681 1.469 0.14177439298
5 2.013 NA NA -0.926 NA NA 3 -50.152 108.150 1.938 0.11211891747
17 1.936 NA NA NA NA 0.971 3 -50.671 109.188 2.976 0.06672651142
29 1.700 NA NA -0.760 0.321 0.721 5 -46.872 109.198 2.986 0.06639978578
23 1.754 NA 0.267 -0.872 NA 0.798 5 -47.065 109.585 3.373 0.05471795199
18 1.797 0.595 NA NA NA 1.037 4 -49.351 110.035 3.823 0.04369206674
13 1.956 NA NA -0.878 0.407 NA 4 -49.662 110.657 4.444 0.03202041698
6 1.964 0.398 NA -0.916 NA NA 4 -49.717 110.767 4.554 0.03031062811
24 1.611 0.611 0.396 -0.873 NA 0.891 6 -45.263 110.927 4.714 0.02797872063
25 1.868 NA NA NA 0.360 0.902 4 -50.087 111.508 5.295 0.02092303537
7 2.004 NA 0.118 -0.952 NA NA 4 -50.094 111.522 5.310 0.02077535264
31 1.641 NA 0.476 -0.884 0.499 0.758 6 -45.830 112.060 5.847 0.01587836355
30 1.656 0.530 NA -0.757 -0.002 0.832 6 -46.113 112.627 6.414 0.01195950591
19 1.938 NA -0.025 NA NA 0.967 4 -50.669 112.671 6.459 0.01169410704
1 2.360 NA NA NA NA NA 2 -54.410 113.677 7.465 0.00707234126
15 1.909 NA 0.370 -0.950 0.619 NA 5 -49.208 113.871 7.658 0.00642054976
20 1.789 0.608 0.071 NA NA 1.046 5 -49.334 114.122 7.910 0.00566196129
26 1.797 0.585 NA NA 0.014 1.034 5 -49.350 114.156 7.943 0.00556800220
8 1.943 0.476 0.221 -0.961 NA NA 5 -49.531 114.517 8.304 0.00464786372
14 1.949 0.219 NA -0.890 0.284 NA 5 -49.574 114.603 8.391 0.00445190818
9 2.284 NA NA NA 0.607 NA 3 -53.815 115.475 9.263 0.00287808609
27 1.853 NA 0.121 NA 0.392 0.912 5 -50.043 115.540 9.328 0.00278593207
2 2.312 0.531 NA NA NA NA 3 -54.042 115.929 9.717 0.00229363542
3 2.370 NA -0.243 NA NA NA 3 -54.333 116.513 10.301 0.00171310123
32 1.600 0.481 0.446 -0.876 0.193 0.856 7 -45.155 116.755 10.542 0.00151813805
16 1.905 0.169 0.355 -0.956 0.514 NA 6 -49.159 118.718 12.505 0.00056887139
10 2.282 0.128 NA NA 0.539 NA 4 -53.801 118.934 12.722 0.00051047033
11 2.288 NA -0.048 NA 0.593 NA 4 -53.812 118.957 12.744 0.00050485639
28 1.788 0.583 0.078 NA 0.035 1.038 6 -49.331 119.062 12.849 0.00047898391
4 2.321 0.493 -0.143 NA NA NA 4 -54.017 119.368 13.156 0.00041097803
12 2.287 0.134 -0.061 NA 0.519 NA 5 -53.796 123.046 16.834 0.00006532418

# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.migra.4 <- get.models(model.set.migra, subset = delta < 4)
avg_model.migra.4 <- model.avg(model.set.migra.4)
summary(avg_model.migra.4)
## 
## Call:
## model.avg(object = model.set.migra.4)
## 
## Component model call: 
## glm.nb(formula = abundance ~ <7 unique rhs>, data = b.migra.hab, 
##      na.action = na.fail, init.theta = <7 unique values>, link = log)
## 
## Component models: 
##     df logLik   AICc delta weight
## 35   4 -47.44 106.21  0.00   0.38
## 135  5 -46.11 107.68  1.47   0.18
## 3    3 -50.15 108.15  1.94   0.14
## 5    3 -50.67 109.19  2.98   0.09
## 345  5 -46.87 109.20  2.99   0.09
## 235  5 -47.07 109.59  3.37   0.07
## 15   4 -49.35 110.04  3.82   0.06
## 
## Term codes: 
## meanturb     s.do    s.sal stnwidth  vegelev 
##        1        2        3        4        5 
## 
## Model-averaged coefficients:  
## (full average) 
##             Estimate Std. Error Adjusted SE z value  Pr(>|z|)    
## (Intercept)  1.79530    0.31154     0.33729   5.323 0.0000001 ***
## s.sal       -0.69562    0.38277     0.39913   1.743    0.0814 .  
## vegelev      0.69379    0.40250     0.42072   1.649    0.0991 .  
## meanturb     0.12919    0.27823     0.28702   0.450    0.6526    
## stnwidth     0.02733    0.12363     0.12978   0.211    0.8332    
## s.do         0.01872    0.10727     0.11393   0.164    0.8695    
##  
## (conditional average) 
##             Estimate Std. Error Adjusted SE z value  Pr(>|z|)    
## (Intercept)   1.7953     0.3115      0.3373   5.323 0.0000001 ***
## s.sal        -0.8102     0.2790      0.3045   2.661    0.0078 ** 
## vegelev       0.8101     0.3081      0.3353   2.416    0.0157 *  
## meanturb      0.5439     0.3168      0.3482   1.562    0.1183    
## stnwidth      0.3214     0.2919      0.3218   0.999    0.3179    
## s.do          0.2672     0.3127      0.3447   0.775    0.4382    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      s.sal vegelev meanturb stnwidth s.do
## Importance:          0.86  0.86    0.24     0.09     0.07
## N containing models:    5     6       2        1        1
migra.ci <- data.frame(confint(avg_model.migra.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4
migra.4.Rsq <- rsquared(model.set.migra.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.migra <- data.frame(avg_model.migra.4$msTable)
avg_model_components4.migra <- cbind(migra.4.Rsq, avg_model_4df.migra)
r = data.frame(Coeff = rownames(avg_model_4df.migra, rep(NA, length(avg_model_components4.migra))))
avg_model_components4.migra <- cbind(avg_model_components4.migra, r)
avg_model_components4.migra <- avg_model_components4.migra[, -c(6, 7)]
# write.csv(avg_model_components4.migra,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachmigratory2.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

Check fit of top model by avg – including all parameters with weight >0.5

avg.migra <- glm.nb(abundance ~ s.sal + vegelev, data = b.migra.hab, na.action = "na.fail")

summary(avg.migra)  # AIC 102.88 
## 
## Call:
## glm.nb(formula = abundance ~ s.sal + vegelev, data = b.migra.hab, 
##     na.action = "na.fail", init.theta = 0.9382988358, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.06918  -0.76237  -0.42767   0.04231   1.89455  
## 
## Coefficients:
##             Estimate Std. Error z value       Pr(>|z|)    
## (Intercept)   1.7763     0.2818   6.303 0.000000000292 ***
## s.sal        -0.7913     0.2645  -2.992        0.00277 ** 
## vegelev       0.7525     0.2910   2.586        0.00971 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9383) family taken to be 1)
## 
##     Null deviance: 36.393  on 16  degrees of freedom
## Residual deviance: 16.677  on 14  degrees of freedom
## AIC: 102.88
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.938 
##           Std. Err.:  0.332 
## 
##  2 x log-likelihood:  -94.879
vif(avg.migra)  #very low
##    s.sal  vegelev 
## 1.008413 1.008413
### Plot residuals vs. fitted values
plot(fitted(avg.migra), resid(avg.migra), main = "Averaged Migratory GLM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = b.migra.hab$abundance, y = resid(avg.migra), main = "Q-Q Averaged Migratory GLM")
qqline(resid(avg.migra), col = 2)

Resident beach seine model

response = average abundance of residents [per beach seine site]
This is 4th of 4 beach models for (1)Chinook, (2)res, (3)other migratory species and (4)resident species

Load all data - standardize variables, create subgroups from beach seine data

### Grab just resident catch - automatically removes any 0s or unidentified
### species.  4
b.4 <- b[which(b$class == "resident"), ]

### standardize continuous variables
b.4$s.temp <- standardize(b.4$Temp.surf, centerFun = mean, scaleFun = sd)
b.4$s.sal <- standardize(b.4$Sal.surf, centerFun = mean, scaleFun = sd)
b.4$s.do <- standardize(b.4$DOmg.surf, centerFun = mean, scaleFun = sd)
b.4$s.pH <- standardize(b.4$pH.surf, centerFun = mean, scaleFun = sd)
b.4$s.J.date <- standardize(b.4$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared
b.4$j2 <- b.4$s.J.date^2

# 1 summarize by site-day
b.res <- ddply(b.4, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal, 
    s.do, s.pH), summarize, abundance = sum(abundance))

# plot abundance over julian day to see distribution:
plot(abundance ~ J.date, data = b.res)

Habitat variables

## Add marsh habitat variables for each site
b.res.hab <- b.res
b.res.hab$stnwidth <- hab2[match(b.res.hab$Site, hab$Site), 2]
b.res.hab$vegelev <- hab2[match(b.res.hab$Site, hab$Site), 3]
b.res.hab$shtdensity <- hab2[match(b.res.hab$Site, hab$Site), 4]
b.res.hab$shtdenhigh <- hab2[match(b.res.hab$Site, hab$Site), 5]
b.res.hab$tcelev <- hab2[match(b.res.hab$Site, hab$Site), 6]
b.res.hab$angbank <- hab2[match(b.res.hab$Site, hab$Site), 7]
b.res.hab$meanturb <- hab2[match(b.res.hab$Site, hab$Site), 8]

## standardize habitat variables using 'standardize' function from robustHD
## package
b.res.hab$stnwidth <- standardize(b.res.hab$stnwidth, centerFun = mean, scaleFun = sd)
b.res.hab$stnwidth <- as.numeric(b.res.hab$stnwidth)
b.res.hab$vegelev <- standardize(b.res.hab$vegelev, centerFun = mean, scaleFun = sd)
b.res.hab$vegelev <- as.numeric(b.res.hab$vegelev)
b.res.hab$shtdensity <- standardize(b.res.hab$shtdensity, centerFun = mean, 
    scaleFun = sd)
b.res.hab$shtdensity <- as.numeric(b.res.hab$shtdensity)
b.res.hab$shtdenhigh <- standardize(b.res.hab$shtdenhigh, centerFun = mean, 
    scaleFun = sd)
b.res.hab$shtdenhigh <- as.numeric(b.res.hab$shtdenhigh)
b.res.hab$tcelev <- standardize(b.res.hab$tcelev, centerFun = mean, scaleFun = sd)
b.res.hab$tcelev <- as.numeric(b.res.hab$tcelev)
b.res.hab$angbank <- standardize(b.res.hab$angbank, centerFun = mean, scaleFun = sd)
b.res.hab$angbank <- as.numeric(b.res.hab$angbank)
b.res.hab$meanturb <- standardize(b.res.hab$meanturb, centerFun = mean, scaleFun = sd)
b.res.hab$meanturb <- as.numeric(b.res.hab$meanturb)

VIF for collinearity of habitat variables

Assess variance inflation factors

# alias function identifies covariates that are multiples of each other - in
# this case have some habitat parameters causing issues.
alias(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + stnwidth + 
    vegelev + shtdenhigh + shtdensity + tcelev + angbank + meanturb, data = b.res.hab))
## Model :
## abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + 
##     stnwidth + vegelev + shtdenhigh + shtdensity + tcelev + angbank + 
##     meanturb
## 
## Complete :
##          (Intercept)         s.J.date            j2                 
## tcelev                     0                   0                   0
## angbank                    0                   0                   0
## meanturb                   0                   0                   0
##          Year                s.temp              s.sal              
## tcelev                     0                   0                   0
## angbank                    0                   0                   0
## meanturb                   0                   0                   0
##          s.do                s.pH                stnwidth           
## tcelev                     0                   0            877/1622
## angbank                    0                   0      -293738/568293
## meanturb                   0                   0          8163/23578
##          vegelev             shtdenhigh          shtdensity         
## tcelev   158272651/348415924   -6656702/13163509            991/3356
## angbank        141665/188257           -948/2101    -955547/14904747
## meanturb      -248443/329042        -55259/84438      -185123/372691

# We determined hab variables for Resident should be reduced to shoot
# density, angbank, and tcelev based on significance and collinearity when
# doing model exploration.
vif(lm(abundance ~ s.J.date + j2 + Year + s.temp + s.sal + s.do + s.pH + shtdensity + 
    tcelev + angbank, data = b.res.hab))  ##all less than 5
##   s.J.date         j2       Year     s.temp      s.sal       s.do 
##   3.835484   2.948591   1.521699   5.571724   2.106828   1.408997 
##       s.pH shtdensity     tcelev    angbank 
##   1.238688   1.143590   1.935296   1.789858


## Pearson Corr with all vars
year <- b.res.hab$Year
jday <- b.res.hab$s.J.date
j2 <- b.res.hab$j2
temp <- b.res.hab$s.temp
do <- b.res.hab$s.do
ph <- b.res.hab$s.pH
sal <- b.res.hab$s.sal
veg <- b.res.hab$vegelev
turb <- b.res.hab$meanturb
width <- b.res.hab$stnwidth
shootden <- b.res.hab$shtdensity
shtdenhigh <- b.res.hab$shtdenhigh
elev <- b.res.hab$tcelev
angbank <- b.res.hab$angbank

habcovar <- cbind.data.frame(year, jday, j2, temp, do, ph, sal, veg, turb, width, 
    shootden, shtdenhigh, elev, angbank)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot for all variables in beach seine resident fish dataset")

ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_residentmarsh.pdf")

Model selection

Full model

Note shoot density came out as a stronger habitat variable than turbidity in this model (delta AIC ~5, stronger weight to top model when run without shtden or meanturb). However, it was strongly correlated to turbidity (0.779) and also highly correlated to stnwidth (0.537). To avoid removing two habitat variables, and to remain consistent with the other global models, we decided not to use shoot density. Note that temperature was very highly correlated to Julian day, J2, and salinity in this data set, and caused several of the VIF values to inflate. Once removed, all VIFs dropped to below 2 and the AIC improved, so it was removed from the global model.

## Resident: full model with top habitat variables
res <- glm.nb(abundance ~ s.J.date + j2 + Year + s.sal + s.do + s.pH + meanturb + 
    vegelev + stnwidth, data = b.res.hab, na.action = "na.fail")
summary(res)  #749.43, VIF below 2  
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.sal + s.do + 
##     s.pH + meanturb + vegelev + stnwidth, data = b.res.hab, na.action = "na.fail", 
##     init.theta = 1.206401257, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6124  -0.8895  -0.3893   0.0975   3.2165  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.98477    0.19314  20.631  < 2e-16 ***
## s.J.date     0.79278    0.11210   7.072 1.53e-12 ***
## j2          -0.48207    0.08914  -5.408 6.37e-08 ***
## Year         1.06309    0.24228   4.388 1.15e-05 ***
## s.sal       -0.07968    0.12752  -0.625 0.532084    
## s.do         0.18154    0.12509   1.451 0.146708    
## s.pH        -0.33039    0.11572  -2.855 0.004301 ** 
## meanturb     0.48323    0.13005   3.716 0.000203 ***
## vegelev      0.03291    0.13490   0.244 0.807294    
## stnwidth    -0.17553    0.13816  -1.270 0.203934    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.2064) family taken to be 1)
## 
##     Null deviance: 196.121  on 77  degrees of freedom
## Residual deviance:  85.743  on 68  degrees of freedom
## AIC: 749.43
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.206 
##           Std. Err.:  0.183 
## 
##  2 x log-likelihood:  -727.426
res2 <- glm.nb(abundance ~ s.J.date + j2 + Year + s.sal + s.do + s.pH + shtdensity + 
    vegelev, data = b.res.hab, na.action = "na.fail")
summary(res2)  #744.11, VIF still below 2
## 
## Call:
## glm.nb(formula = abundance ~ s.J.date + j2 + Year + s.sal + s.do + 
##     s.pH + shtdensity + vegelev, data = b.res.hab, na.action = "na.fail", 
##     init.theta = 1.244095276, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7463  -1.0202  -0.3433   0.3050   3.1930  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.96891    0.18948  20.946  < 2e-16 ***
## s.J.date     0.75583    0.11092   6.814 9.49e-12 ***
## j2          -0.47398    0.08801  -5.385 7.23e-08 ***
## Year         1.00505    0.23662   4.247 2.16e-05 ***
## s.sal       -0.06401    0.12392  -0.517   0.6055    
## s.do         0.22313    0.12327   1.810   0.0703 .  
## s.pH        -0.33959    0.11443  -2.968   0.0030 ** 
## shtdensity  -0.48810    0.11047  -4.418 9.94e-06 ***
## vegelev     -0.21295    0.11279  -1.888   0.0590 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.2441) family taken to be 1)
## 
##     Null deviance: 201.965  on 77  degrees of freedom
## Residual deviance:  84.981  on 69  degrees of freedom
## AIC: 744.21
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.244 
##           Std. Err.:  0.190 
## 
##  2 x log-likelihood:  -724.206
vif(res)
## s.J.date       j2     Year    s.sal     s.do     s.pH meanturb  vegelev 
## 1.366545 1.911245 1.280196 1.986694 1.424802 1.269219 1.457012 1.601792 
## stnwidth 
## 1.675253
# Check model assumptions - note that several parameters have non ideal
# relationships. Also note that there are several outliers, but these are
# real so we do not remove them (despite what this test says).
sjp.glm(res, type = "ma")
## `sjp.glm()` will become deprecated in the future. Please use `plot_model()` instead.
## 1 outliers removed in updated model.
## # A tibble: 2 x 2
##   models     aic
##   <chr>    <dbl>
## 1 original  749.
## 2 updated   718.

## --------------------
## Check significance of terms when they entered the model...
## Anova:
## Warning in anova.negbin(logreg, test = "Chisq"): tests made without re-
## estimating 'theta'

## Analysis of Deviance Table
## 
## Model: Negative Binomial(1.2064), link: log
## 
## Response: abundance
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        77    196.121              
## s.J.date  1   10.356        76    185.765 0.0012906 ** 
## j2        1   56.016        75    129.749 7.187e-14 ***
## Year      1   16.037        74    113.711 6.210e-05 ***
## s.sal     1    0.001        73    113.711 0.9807157    
## s.do      1    1.175        72    112.536 0.2783492    
## s.pH      1   10.586        71    101.950 0.0011394 ** 
## meanturb  1   14.465        70     87.485 0.0001428 ***
## vegelev   1    0.237        69     87.248 0.6261642    
## stnwidth  1    1.505        68     85.743 0.2199672    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# See pseudo R2 for full model
rsquared(res, aicc = TRUE)
##    Class            Family Link  n R.squared
## 1 negbin Negative Binomial  log 78  0.562805

Model selection - dredge

Then continued to dredge function to determine optimal model using AICc (res has 17 observations).

# Generate model set
model.set.res <- dredge(res)
# Create model selection table
model_table.res <- model.sel(model.set.res)
options(scipen = 7)
names(model_table.res) <- c("(Int)", "Julian day^2", "Mean turb", "DO", "Julian day", 
    "pH", "Sal", "Channel width", "Veg. elev.", "Year", "df", "LL", "AICc", 
    "delta", "weight")
kable(head(model_table.res, n = 100), digits = 3)
(Int) Julian day^2 Mean turb DO Julian day pH Sal Channel width Veg. elev. Year df LL AICc delta weight
284 4.129 -0.529 0.442 NA 0.740 -0.301 NA NA NA 0.920 7 -365.573 746.746 0.000 0.182920560692
288 4.080 -0.520 0.418 0.168 0.778 -0.349 NA NA NA 0.954 8 -364.827 747.741 0.994 0.111266844172
348 4.087 -0.529 0.484 NA 0.756 -0.285 NA -0.149 NA 1.000 8 -364.841 747.768 1.022 0.109746806824
352 4.024 -0.516 0.461 0.188 0.798 -0.335 NA -0.165 NA 1.046 9 -363.912 748.471 1.725 0.077230890598
316 4.074 -0.491 0.452 NA 0.736 -0.301 -0.087 NA NA 0.959 8 -365.339 748.765 2.019 0.066654400522
412 4.118 -0.530 0.430 NA 0.742 -0.303 NA NA -0.043 0.945 8 -365.501 749.089 2.343 0.056694729815
320 4.023 -0.480 0.429 0.170 0.774 -0.349 -0.090 NA NA 0.991 9 -364.575 749.796 3.050 0.039808074651
416 4.055 -0.519 0.394 0.186 0.784 -0.357 NA NA -0.074 0.997 9 -364.621 749.889 3.143 0.038007892528
380 4.043 -0.497 0.491 NA 0.751 -0.286 -0.075 -0.144 NA 1.029 9 -364.659 749.966 3.219 0.036577090674
476 4.091 -0.528 0.502 NA 0.757 -0.280 NA -0.172 0.042 0.990 9 -364.793 750.233 3.487 0.031996820852
268 4.085 -0.488 0.465 NA 0.734 NA NA NA NA 1.027 6 -368.526 750.235 3.489 0.031968121608
332 4.021 -0.488 0.517 NA 0.751 NA NA -0.198 NA 1.142 7 -367.362 750.324 3.578 0.030577406544
384 3.984 -0.485 0.468 0.187 0.793 -0.335 -0.073 -0.159 NA 1.070 10 -363.740 750.763 4.017 0.024549895481
480 4.026 -0.516 0.466 0.186 0.798 -0.333 NA -0.172 0.012 1.043 10 -363.908 751.100 4.354 0.020740718087
444 4.072 -0.495 0.444 NA 0.737 -0.303 -0.080 NA -0.027 0.972 9 -365.313 751.273 4.527 0.019020837136
448 4.014 -0.486 0.408 0.183 0.778 -0.356 -0.074 NA -0.058 1.019 10 -364.459 752.201 5.454 0.011964306605
336 3.976 -0.475 0.509 0.104 0.774 NA NA -0.214 NA 1.186 8 -367.096 752.278 5.532 0.011509420075
300 4.030 -0.451 0.477 NA 0.732 NA -0.080 NA NA 1.068 7 -368.347 752.294 5.548 0.011415914866
508 4.042 -0.490 0.520 NA 0.752 -0.278 -0.089 -0.179 0.065 1.017 10 -364.554 752.391 5.645 0.010877766864
272 4.063 -0.481 0.458 0.064 0.748 NA NA NA NA 1.048 7 -368.425 752.449 5.703 0.010564950833
460 4.033 -0.488 0.549 NA 0.755 NA NA -0.239 0.080 1.111 8 -367.202 752.490 5.744 0.010350615802
364 3.980 -0.459 0.525 NA 0.748 NA -0.065 -0.195 NA 1.172 8 -367.233 752.553 5.806 0.010033063566
396 4.073 -0.488 0.456 NA 0.734 NA NA NA -0.034 1.052 7 -368.485 752.571 5.824 0.009942402590
512 3.985 -0.482 0.483 0.182 0.793 -0.330 -0.080 -0.176 0.033 1.063 11 -363.713 753.426 6.680 0.006482096312
492 3.980 -0.445 0.572 NA 0.753 NA -0.094 -0.248 0.111 1.138 9 -366.955 754.557 7.810 0.003683410824
304 4.007 -0.443 0.469 0.064 0.746 NA -0.081 NA NA 1.088 8 -368.243 754.572 7.826 0.003654809769
368 3.938 -0.447 0.517 0.102 0.771 NA -0.063 -0.210 NA 1.212 9 -366.976 754.599 7.853 0.003606288205
464 3.991 -0.476 0.536 0.094 0.775 NA NA -0.246 0.066 1.157 9 -366.992 754.630 7.884 0.003550787903
428 4.029 -0.453 0.472 NA 0.732 NA -0.075 NA -0.015 1.077 8 -368.340 754.768 8.021 0.003314812243
400 4.040 -0.479 0.443 0.077 0.750 NA NA NA -0.050 1.088 8 -368.345 754.777 8.030 0.003299759334
496 3.945 -0.437 0.559 0.087 0.772 NA -0.089 -0.254 0.096 1.177 10 -366.775 756.834 10.087 0.001179934982
432 4.001 -0.447 0.459 0.072 0.747 NA -0.071 NA -0.030 1.108 9 -368.217 757.080 10.334 0.001042883619
414 4.141 -0.542 NA 0.294 0.761 -0.419 NA NA -0.193 0.984 8 -369.676 757.439 10.693 0.000871613823
286 4.235 -0.554 NA 0.262 0.744 -0.412 NA NA NA 0.856 7 -371.167 757.934 11.188 0.000680572339
282 4.317 -0.572 NA NA 0.696 -0.338 NA NA NA 0.820 6 -372.714 758.611 11.864 0.000485219072
410 4.250 -0.568 NA NA 0.709 -0.343 NA NA -0.163 0.928 7 -371.651 758.902 12.155 0.000419546007
28 4.642 -0.555 0.396 NA 0.631 -0.371 NA NA NA NA 6 -372.956 759.094 12.348 0.000381018984
478 4.144 -0.542 NA 0.289 0.759 -0.423 NA 0.043 -0.209 0.973 9 -369.623 759.893 13.147 0.000255537591
446 4.145 -0.546 NA 0.294 0.762 -0.419 0.008 NA -0.194 0.982 9 -369.674 759.996 13.249 0.000242779941
350 4.221 -0.553 NA 0.272 0.750 -0.409 NA -0.050 NA 0.883 8 -371.074 760.235 13.489 0.000215355967
318 4.217 -0.539 NA 0.263 0.744 -0.413 -0.033 NA NA 0.866 8 -371.138 760.362 13.616 0.000202111334
314 4.300 -0.559 NA NA 0.696 -0.338 -0.030 NA NA 0.832 7 -372.689 760.978 14.232 0.000148554064
32 4.622 -0.549 0.373 0.097 0.652 -0.404 NA NA NA NA 7 -372.695 760.991 14.245 0.000147604798
346 4.313 -0.573 NA NA 0.699 -0.336 NA -0.022 NA 0.832 7 -372.696 760.993 14.246 0.000147479176
474 4.252 -0.566 NA NA 0.706 -0.349 NA 0.063 -0.186 0.912 8 -371.537 761.162 14.415 0.000135535662
156 4.637 -0.552 0.415 NA 0.634 -0.357 NA NA 0.066 NA 7 -372.809 761.217 14.471 0.000131818091
442 4.250 -0.568 NA NA 0.709 -0.343 -0.001 NA -0.163 0.929 8 -371.651 761.389 14.642 0.000120992507
60 4.654 -0.568 0.394 NA 0.635 -0.370 0.028 NA NA NA 7 -372.936 761.472 14.725 0.000116061701
92 4.642 -0.555 0.401 NA 0.632 -0.372 NA -0.016 NA NA 7 -372.948 761.496 14.749 0.000114689361
266 4.314 -0.536 NA NA 0.680 NA NA NA NA 0.869 5 -375.831 762.495 15.748 0.000069592784
510 4.147 -0.545 NA 0.289 0.759 -0.423 0.006 0.042 -0.209 0.971 10 -369.622 762.528 15.781 0.000068455690
382 4.207 -0.542 NA 0.272 0.749 -0.410 -0.027 -0.047 NA 0.890 9 -371.055 762.756 16.010 0.000061064849
394 4.226 -0.531 NA NA 0.685 NA NA NA -0.161 1.016 6 -374.899 762.980 16.234 0.000054592919
160 4.622 -0.548 0.391 0.086 0.652 -0.390 NA NA 0.050 NA 8 -372.615 763.316 16.570 0.000046148090
378 4.297 -0.560 NA NA 0.698 -0.337 -0.029 -0.020 NA 0.842 8 -372.674 763.435 16.688 0.000043493776
96 4.622 -0.550 0.379 0.099 0.654 -0.406 NA -0.024 NA NA 8 -372.678 763.444 16.697 0.000043298394
220 4.638 -0.553 0.447 NA 0.640 -0.353 NA -0.076 0.108 NA 8 -372.682 763.452 16.705 0.000043123246
64 4.632 -0.560 0.372 0.096 0.656 -0.402 0.022 NA NA NA 8 -372.683 763.452 16.706 0.000043112910
188 4.640 -0.555 0.414 NA 0.635 -0.357 0.006 NA 0.065 NA 8 -372.808 763.703 16.956 0.000038043331
506 4.252 -0.566 NA NA 0.706 -0.349 -0.001 0.063 -0.186 0.912 9 -371.537 763.722 16.975 0.000037683780
124 4.655 -0.569 0.399 NA 0.637 -0.370 0.029 -0.018 NA NA 8 -372.926 763.939 17.193 0.000033795603
270 4.264 -0.520 NA 0.145 0.705 NA NA NA NA 0.901 6 -375.383 763.948 17.202 0.000033644199
398 4.142 -0.506 NA 0.187 0.715 NA NA NA -0.187 1.083 7 -374.176 763.952 17.206 0.000033580163
330 4.298 -0.538 NA NA 0.685 NA NA -0.059 NA 0.907 6 -375.716 764.614 17.868 0.000024117160
298 4.312 -0.534 NA NA 0.680 NA -0.004 NA NA 0.871 6 -375.830 764.843 18.097 0.000021504137
426 4.246 -0.547 NA NA 0.686 NA 0.036 NA -0.167 1.003 7 -374.865 765.330 18.584 0.000016860157
458 4.226 -0.530 NA NA 0.684 NA NA 0.009 -0.164 1.014 7 -374.896 765.393 18.646 0.000016340838
224 4.622 -0.549 0.422 0.084 0.658 -0.384 NA -0.073 0.090 NA 9 -372.499 765.645 18.898 0.000014406394
12 4.694 -0.516 0.382 NA 0.658 NA NA NA NA NA 5 -377.426 765.686 18.940 0.000014109604
192 4.624 -0.551 0.390 0.086 0.653 -0.390 0.006 NA 0.049 NA 9 -372.614 765.875 19.128 0.000012841836
334 4.231 -0.518 NA 0.167 0.714 NA NA -0.085 NA 0.960 7 -375.152 765.905 19.158 0.000012650753
128 4.633 -0.561 0.378 0.098 0.658 -0.404 0.024 -0.025 NA NA 9 -372.663 765.974 19.227 0.000012221557
252 4.637 -0.552 0.448 NA 0.640 -0.353 -0.003 -0.077 0.108 NA 9 -372.682 766.012 19.265 0.000011991232
430 4.167 -0.529 NA 0.191 0.717 NA 0.051 NA -0.197 1.069 8 -374.112 766.312 19.565 0.000010321217
302 4.263 -0.519 NA 0.145 0.705 NA -0.002 NA NA 0.902 7 -375.382 766.365 19.619 0.000010049514
462 4.140 -0.506 NA 0.189 0.716 NA NA -0.013 -0.182 1.087 8 -374.171 766.430 19.683 0.000009729040
140 4.680 -0.512 0.426 NA 0.668 NA NA NA 0.143 NA 6 -376.761 766.705 19.959 0.000008477039
26 4.767 -0.593 NA NA 0.596 -0.366 NA NA NA NA 5 -377.950 766.734 19.987 0.000008357446
362 4.297 -0.537 NA NA 0.685 NA -0.001 -0.059 NA 0.908 7 -375.715 767.031 20.285 0.000007203186
30 4.712 -0.576 NA 0.199 0.633 -0.433 NA NA NA NA 6 -376.969 767.121 20.375 0.000006884766
44 4.733 -0.558 0.374 NA 0.668 NA 0.084 NA NA NA 6 -377.256 767.695 20.948 0.000005168447
490 4.247 -0.547 NA NA 0.685 NA 0.037 0.011 -0.171 1.000 8 -374.862 767.811 21.065 0.000004877129
16 4.702 -0.521 0.391 -0.038 0.649 NA NA NA NA NA 6 -377.388 767.959 21.213 0.000004528053
76 4.694 -0.516 0.383 NA 0.658 NA NA -0.005 NA NA 6 -377.426 768.035 21.288 0.000004360744
256 4.622 -0.548 0.422 0.084 0.658 -0.384 -0.001 -0.073 0.091 NA 10 -372.499 768.281 21.535 0.000003855348
366 4.234 -0.520 NA 0.167 0.714 NA 0.004 -0.085 NA 0.958 8 -375.152 768.391 21.644 0.000003650020
204 4.679 -0.513 0.478 NA 0.679 NA NA -0.127 0.210 NA 7 -376.442 768.485 21.738 0.000003482261
90 4.757 -0.588 NA NA 0.593 -0.364 NA 0.072 NA NA 6 -377.773 768.728 21.982 0.000003082673
154 4.765 -0.594 NA NA 0.595 -0.378 NA NA -0.055 NA 6 -377.839 768.861 22.115 0.000002885089
494 4.165 -0.528 NA 0.193 0.718 NA 0.050 -0.013 -0.192 1.073 9 -374.108 768.863 22.117 0.000002881381
58 4.794 -0.622 NA NA 0.603 -0.362 0.062 NA NA NA 6 -377.859 768.901 22.155 0.000002827811
144 4.690 -0.518 0.444 -0.063 0.652 NA NA NA 0.150 NA 7 -376.659 768.918 22.172 0.000002803762
158 4.703 -0.575 NA 0.209 0.633 -0.452 NA NA -0.078 NA 7 -376.750 769.100 22.354 0.000002559527
172 4.689 -0.521 0.423 NA 0.669 NA 0.018 NA 0.136 NA 7 -376.755 769.110 22.364 0.000002547133
94 4.709 -0.575 NA 0.190 0.629 -0.429 NA 0.049 NA NA 7 -376.885 769.370 22.624 0.000002236313
62 4.734 -0.600 NA 0.196 0.638 -0.429 0.050 NA NA NA 7 -376.911 769.421 22.675 0.000002180297
299 3.558 NA 0.599 NA 0.669 NA -0.439 NA NA 1.192 6 -378.375 769.932 23.186 0.000001688567
48 4.740 -0.563 0.383 -0.039 0.658 NA 0.084 NA NA NA 7 -377.217 770.033 23.287 0.000001605295
108 4.734 -0.559 0.377 NA 0.668 NA 0.086 -0.014 NA NA 7 -377.251 770.102 23.355 0.000001551472
80 4.702 -0.521 0.392 -0.038 0.649 NA NA -0.002 NA NA 7 -377.388 770.376 23.630 0.000001352465

# Model averaging Version 1: use all models with delta AIC score less than 4
model.set.res.4 <- get.models(model.set.res, subset = delta < 4)
avg_model.res.4 <- model.avg(model.set.res.4)
summary(avg_model.res.4)
## 
## Call:
## model.avg(object = model.set.res.4)
## 
## Component model call: 
## glm.nb(formula = abundance ~ <12 unique rhs>, data = b.res.hab, 
##      na.action = na.fail, init.theta = <12 unique values>, link = log)
## 
## Component models: 
##         df  logLik   AICc delta weight
## 12459    7 -365.57 746.75  0.00   0.22
## 123459   8 -364.83 747.74  0.99   0.14
## 124579   8 -364.84 747.77  1.02   0.13
## 1234579  9 -363.91 748.47  1.72   0.09
## 124569   8 -365.34 748.77  2.02   0.08
## 124589   8 -365.50 749.09  2.34   0.07
## 1234569  9 -364.57 749.80  3.05   0.05
## 1234589  9 -364.62 749.89  3.14   0.05
## 1245679  9 -364.66 749.97  3.22   0.04
## 1245789  9 -364.79 750.23  3.49   0.04
## 1249     6 -368.53 750.23  3.49   0.04
## 12479    7 -367.36 750.32  3.58   0.04
## 
## Term codes: 
##       j2 meanturb     s.do s.J.date     s.pH    s.sal stnwidth  vegelev 
##        1        2        3        4        5        6        7        8 
##     Year 
##        9 
## 
## Model-averaged coefficients:  
## (full average) 
##              Estimate Std. Error Adjusted SE z value  Pr(>|z|)    
## (Intercept)  4.081422   0.180832    0.183840  22.201   < 2e-16 ***
## j2          -0.516061   0.073616    0.074828   6.897   < 2e-16 ***
## meanturb     0.451576   0.117970    0.119888   3.767  0.000165 ***
## s.J.date     0.757814   0.107372    0.109167   6.942   < 2e-16 ***
## s.pH        -0.288963   0.138021    0.139470   2.072  0.038277 *  
## Year         0.979812   0.237464    0.241384   4.059 0.0000493 ***
## s.do         0.057912   0.109888    0.110713   0.523  0.600919    
## stnwidth    -0.056596   0.105597    0.106453   0.532  0.594969    
## s.sal       -0.014969   0.061265    0.062036   0.241  0.809324    
## vegelev     -0.004848   0.050826    0.051565   0.094  0.925098    
##  
## (conditional average) 
##             Estimate Std. Error Adjusted SE z value  Pr(>|z|)    
## (Intercept)  4.08142    0.18083     0.18384  22.201   < 2e-16 ***
## j2          -0.51606    0.07362     0.07483   6.897   < 2e-16 ***
## meanturb     0.45158    0.11797     0.11989   3.767  0.000165 ***
## s.J.date     0.75781    0.10737     0.10917   6.942   < 2e-16 ***
## s.pH        -0.31303    0.11446     0.11635   2.690  0.007137 ** 
## Year         0.97981    0.23746     0.24138   4.059 0.0000493 ***
## s.do         0.17689    0.12585     0.12804   1.382  0.167114    
## stnwidth    -0.16090    0.12214     0.12424   1.295  0.195287    
## s.sal       -0.08513    0.12399     0.12615   0.675  0.499785    
## vegelev     -0.03112    0.12557     0.12749   0.244  0.807126    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Relative variable importance: 
##                      j2   meanturb s.J.date Year s.pH stnwidth s.do s.sal
## Importance:          1.00 1.00     1.00     1.00 0.92 0.35     0.33 0.18 
## N containing models:   12   12       12       12   10    5        4    3 
##                      vegelev
## Importance:          0.16   
## N containing models:    3
res.ci <- data.frame(confint(avg_model.res.4, full = TRUE))

# Get pseudo R squared values for models up to delta < 4
res.4.Rsq <- rsquared(model.set.res.4, aicc = TRUE)

## write tables to .csv for easy comparison and plugging into appendix table
avg_model_4df.res <- data.frame(avg_model.res.4$msTable)
avg_model_components4.res <- cbind(res.4.Rsq, avg_model_4df.res)
r = data.frame(Coeff = rownames(avg_model_4df.res, rep(NA, length(avg_model_components4.res))))
avg_model_components4.res <- cbind(avg_model_components4.res, r)
avg_model_components4.res <- avg_model_components4.res[, -c(6, 7)]
# write.csv(avg_model_components4.res,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_beachresident.csv')

Parameter Plot

The results of model averaging including all top ranked models up to delta AICc 4

Check fit of top model by avg – including all parameters with weight >0.5

avg.res <- glm.nb(abundance ~ Year + s.J.date + j2 + meanturb + s.pH, data = b.res.hab, 
    na.action = "na.fail")

summary(avg.res)
## 
## Call:
## glm.nb(formula = abundance ~ Year + s.J.date + j2 + meanturb + 
##     s.pH, data = b.res.hab, na.action = "na.fail", init.theta = 1.152490818, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6782  -0.9139  -0.3890   0.1030   3.2421  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.12855    0.17335  23.816  < 2e-16 ***
## Year         0.92038    0.22756   4.045 5.24e-05 ***
## s.J.date     0.74024    0.10212   7.249 4.21e-13 ***
## j2          -0.52937    0.06791  -7.795 6.43e-15 ***
## meanturb     0.44222    0.11029   4.010 6.08e-05 ***
## s.pH        -0.30052    0.10907  -2.755  0.00586 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.1525) family taken to be 1)
## 
##     Null deviance: 187.737  on 77  degrees of freedom
## Residual deviance:  85.784  on 72  degrees of freedom
## AIC: 745.15
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.152 
##           Std. Err.:  0.174 
## 
##  2 x log-likelihood:  -731.146
vif(avg.res)
##     Year s.J.date       j2 meanturb     s.pH 
## 1.086656 1.103259 1.089915 1.007426 1.082319
### Plot residuals vs. fitted values
plot(fitted(avg.res), resid(avg.res), main = "Averaged Resident GLM", xlab = "Fitted", 
    ylab = "Pearson residuals")

## q plot
qqnorm(x = b.res.hab$abundance, y = resid(avg.res), main = "Q-Q Averaged Resident GLM")
qqline(resid(avg.res), col = 2)